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Preface 

This  report  is  a  summary  of  my  study  of  State  Space  optimization 
techniques.  My  purpose  was  to  demonstrate  the  use  of  these  techniques 
by  applying  them  to  a  typical  adaptive  control  system. 

This  study  was  suggested  by  Capt.  F.  M.  Brown,  an  Assistant 
Professor  in  the  Department  of  Electrical  Engineering  at  tbe  Air 
Force  Institute  of  Technology.  The  study  was  further  prompted  by 
the  growing  use  of  State  Space  techniques  in  the  field  of  control 
systems  analysis. 

I  would  like  to  acknowledge  the  support  and  assistance  given 
me  by  Capt.  R.  A.  Hannen,  my  thesis  advisor  and  Assistant  Professor 
in  the  Department  of  Electrical  Engineering  at  the  Institute  of 
Technology. 

Greatful  acknowledgement  is  also  give  to  Mr.  P.  B.  Brentani 
of  the  Minneapolis-Honeywell  Regulator  Company  for  the  helful  ex¬ 
planations  and  information  which  he  sent  me. 

Finally,  I  would  like  to  point  out  that  the  mathematical  derivations 
presented  in  Chapters  III  and  IV  assume  a  knowledge  of  matrix  algebra 
on  the  part  of  the  reader. 

Willian)  B.  Goggins,  Jr. 
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Abstract 


Interest  in  the  field  of  time  optimization  of  control  systems 
has  led  to  the  development  of  a  great  number  of  theories  concerning 
this  subject.  Many  of  these  have  been  .found  difficult  to  apply  to 
a>  practical  system.  However,  Ho  and  Brentani  have  developed  a  way 
to  compute  the  time  optimum  driving  function  to  be  applied  to  a 
linear  controller  by  the  method  of  Steepest  Descent, 

The  method  of  Ho  and  Brentani  was  applied  to  the  linear  con¬ 
troller  of  a  practical  relay  servo  loop  in  order  to  compute  the  time 
optimum  driving  function  for  this  loop.  The  relay  servo  loop  was 
also  simulated  on  an  analogue  computer^  the  results  of  simulation 
were  compared  with  the  time  optimum  driving  function*  It  was  found 
that^the  two  results  correlated  well* 
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I.  Introduction 

One  of  the  most  interesting  and  challenging  problems  in  control 
system  engineering  is  the  application  of  optimization  techniques 
to  an  automatic  control  system.  Throughout  the  years  many  optimization 
techniques  have  been  employed  and  optimizations  of  several  different 
criteria  have  been  attempted.  This  report  will  be  concerned  with 
time  optimization  of  a  system  using  State  Space  techniques.  By  time 
optimization  is  meant  that  a  system  is  constructed  or  operated  in  such 
a  way  as  to  drive  any  system  error  to  zero  in  the  least  time.  The  system 
to  which  these  optimization  techniques  will  be  applied  is  an  adaptive 
control  system  developed  by  the  Minneapolis-Honeywell  Regulator 
Company.  This  system  is  set  forth  in  WADC  Technical  Report  57-349  (Ref  BsB5). 
Since  the  concepts  of  adaptive  control,  of  time  optimization  and  of 
State  Space  techniques  are  not  easily  definable  in  a  sentence  or  two 
a  further  general  discussion  of  the  three  will  be  one  of  the  subjects 
'of  this  chapter.  However,  a  brief  statement  of  the  problem  is  in  order 
at  this  point. 

Statement  of  the  Problem  and  Purpose 

The  basic  problem  considered  in  this  study  is  the  determination 
of  a  time  optimum  driving  function  to  be  applied  to  the  linear  con¬ 
troller  of  an  aircraft.  A  comparison  of  this  optimum  driving  function 
to  that  actually  used  during  control  of  the  same  aircraft  will  then 
be  made.  From  this  comparison  conclusions  will  be  drawn 'pertaihing-to 
the  validity  of  the  optimum  solution  and  to  the  effectiveness  of  the 
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physical  system.  With  this  in  mind,  it  might  be  well  to  consider, 
very  generally,  the  concepts  of  adaptive  control,  of  time  optimi¬ 
zation,  and  of  State. 

The  Adaptive  Concept 

First,  a  brief  discussion  of  the  adaptive  concept  is  in  order. 

In  a  WADE  Technical  Report  an  adaptive  control  system  is  defined  as 
one  "which  has  the  capability  of  changing  its  parameters  through  an 
internal  process  of  measurement,  evaluation,  and  adjustment,  to  adapt 
to  a  changing  environment,  either  external  or  internal,  to  the  vehicle 
under  control  (Ref  2  )•" 

One  way  in  which  such  an  adaptive  system  may  be  realized  is, 
shown  in  Figure  1. 


Model 

6L 

G  =  1 

,  ^0 

F igure  1 

A  Hypothetical  Adaptive  Control  System 

The  command  signal  ©c  is  fed  into  the  model,  the  response  of  which 
is  the  response  desired  for  the  overall  system.  The  rest  of  the  system 
is  then  designed  to  have  a  transfer  function  of  unity.  For  example. 
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an  aircraft  adaptive  autopilot  would  consist  of  the  following  two  basic 
elements;  (l)  a  model,  which  would  have  the  desired  transfer  function 
of  the  autopilot-aircraft  combination,  and  (2)  the  actual  autopilot- 
aircraft  combination  which  would  be  designed  to  have  a  transfer  function 
of  unity  at  all  times.  Thus  the  output  of  the  entire  system  is  that 
of  the  model,  and  the  desired  response  is  obtained.  However, 
since  G  includes  the  aircraft,  the  dynamics  of  which  are  constantly 
changing  throughout  its  flight  envelope,  other  parameters  of  G  must 
be  constantly  changed  in  order  to  keep  G  as  close  to  unity  as  possible. 

The  changing  of  these  parameters  through  internal  measurement,  evaluation, 
and  adjustment  is  the  adaptive  process. 

Time  Optimization  and  the  Adaptive  System 

Obviously  it  is  impossible  to  physically  realize  a  G  of  unity. 
However,  if  a  system  could  be  constructed  which  would  respond  in  minimum 
time  to  the  output  of  a  model,  a  G  approximating  unity  would  be  realized. 

A  general  theory  covering  this  problem  was  developed  in  a  series  of 
papers  by  three  Russians,  Pontriagin,  Boltyanskii,  and  Gramkrelidze, 
and  is  known  as  the  maximum  principal  of  Pontriagin,  (Ref  11:863) 
Application  of  the  maximum  principal  to  a  linear  system  shows  that 
if  the  driving  function  of  the  system  is  "bang-bang",  that  is,  at 
all  times  as  large,  in  one  direction  or  the  other,  as  the  stops 
permit,  then  the  system  will  be  time  optimum.  Such  a  driving  function 
can  be  obtained  from  a  relay.  A  further  obvious  limitation  is  that 
a  proper  switching  sequence  for  the  relay  must  be  provided. 
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An  adaptive  system  based  on  the  maximum  principal  might  therefore 
be  constructed  as  in  Figure  2. 


Model 

■■ 

Switching 

Relay 

Servo 

and 

1 

1 

Logic 

Aircraft 

A 

F igure  2 

Practical  Adaptive  System 

In  this  case  the  output  of  the  model  is  fed  into  a  loop  which  consists 
of  a  switching  logic,  a  relay,  and  a  servo  and  aircraft,  all  enclosed 
by  uni^y  feedback.  The  portion  of  the 'system  of  Figure  2  enclosed  by 
the  feedback  loop  will  be  termed  the  relay  servo  loop  and  will  be 
referred  to  as  such  throughout  the  remainder  of  this  report. 

The  question  of  when  to  switch  the  relay  to  obtain  the  optimum 
response  for  the  above  system  must  still  be  answered.  To  date,  a 
practical  optimum  switching  logic  has  not  been  synthesized.  However, 
for  the  linear  portion  of  the  system  consisting  of  the  servo  and 
aircraft  one  can  compute  relay  switching  times  to  drive  a  given  error 
to  zero  in  minimum  time.  This  is  done  by  a  method  of  successive 
approximations  which  will  be  treated  in  detail  in  Chapters  III  and  IV» 
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The  Concept  of  State 

A  third  consideration  is  the  concept  of  State,  since  techniques 
associated  with  this  concept  will  be  used  to  determine  the  optimum 
switching  times  for  the  aforementioned  relay  (Ref  11:859-860) » 

Very  generally,  one  might  consider  that  the  variables  of  a  system 
consist  of  three  types  of  time  varying  function; 

1)  a  controllable  (input)  variable which  can  be  changed  at  will 
during  the  control  interval 

2)  an  initially  controllable  (state)  variable  TCft) which  can  be  chosen 

at  but  thereafter  will  act  as  a  function  of  M*),  the  input; 

3)  an  observable  (output)  variable  y(t). 

When  the  first  two  of  these  three  functions  are  known  the  system  can 
be  characterized  as  a  function  of  both  of  them.  Also,  the  output  can 
be  uniquely  determined  as  a  function  of  both.  This  concept  may  be 
expressed  in  equation  form: 


I 

yw  ^  (2) 


In  the  case  of  a  system  of  differential  equations  they  assume  the  form 


'X(^)  -  I  CxC^)j 
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A  particular  case  of  the  above  is  the  matrix  equations 
i  ^  X  +  5"^ 


5  (6) 

where  ^ ,  Q ,  and  0  are  non-singular  constant  matrices* 

Any  linear  transfer  function  which  can  be  expressed  in  differential 
equation  form  can  be  written  in  the  form  of  (5)  by  a  suitable  trans¬ 
formation.  The  resulting  equation  is'  then  the  state  equation  of  the 
system. 

A  specific  example  of  a  tranformation  from  a  transfer  function 
to  State  Space  notation  is  in  order.  Consider  the  transfer  function 


5*4  oTj  +  b 


The  associated  differential  equation  is 


(7) 


cV-e, 


-+ 


h  Co 


Let 


(8) 


~  (10) 


^  z.  7^, 
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Substituting  (9),  (10),  and  (11)  into  (8)  yields 

Xt-  -  ~  a.  -h  (12) 

In  matrix  notation  the  equations  are  expressed  as 


Several  points  of  interest  of  the  State  Space  notation  should 
now  be  explained.  If,  at  any  instant  "ic,  and  are  plotted  on  mutually 
perpendicular  axes,  the  resultant  is  a  two  dimensional  vector.  This 
can  be  extended  of  course  for  an  N  dimensional  system  into  N  dimensional 
Euclidian  space;  this  fact  is  the  origin  of  the  terms  State  Space  and 
State  Vector, 

As  time  is  allowed  to  vary,  the  State  Vector  then  describes  a  path 
in  the  N  dimensional  space.  Noting  that  the  components  of  the  State 
Vector  are  the  succesive  derivatives  of  the  first  component,  one  sees 
that  for  a  second  order  system  this  plot  in  a  plane  is  the  familiar 
phase  plane  which  is  quite  popular  for  the  analysis  of  non-linear 
systems.  It  seems  reasonable  that  the  phase  plane  analysis  could 
somehow  be  extended  to  N  dimensional  phase  space. 

However,  State  Space  analysis  is  useful  for  other  reasons. 

One  of  the  more  important  of  these  of  the.  transformation , of ■  the 
system  equations  into  matrix  form  which  facilitates  a  further  trans¬ 
formation  from  one  coordinate  system  to  another,  that  is  from  N 
coordinate  space  to  M  coordinate  space.  The  use  of  matrices  also 
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facilitates  programming  of  a  digital  computer.  In  view  of  the  above 
discussion  a  detailed  statement  of  the  scope  of  investigation  and 
plan  of  developement  can  be  made* 

Scope  of  Investigation  and  Plan  of  Development 

.  An  adaptive  control  system  employing  a  model  followed  by  a  relay 
servo  loop  was  designed,  built  and  flight  tested  in  an  r94C  aircraft 
by  the  Minneapolis-Honeywell  Regulator  Company  in  1957  as  mentioned 
on  page  1  of  this  report.  The  system  was.  a  pitch  axis  controller.  This 
report  will  analyze  the  system  by  State  Space  techniques  in  order  to 
compute  the  time  optimum  driving  function  for  the  linear  portion  of 
the  relay  servo  loop.  A  development  of  these  techniques  will’ first 
be  presented,  however.  The  results  of  an  analogue  computer  simulation 
of  the  relay  servo  loop  will  then  be  used  as  a  basis  for  comparison. 

In  Chapter  II  a  brief  qualitative  explanation  of  the  system  will 
be  presented  in  order  to  explain  the  operation  of  the  basic  system 
and  to  set  the  stage  for  the  analysis  of  the  succeeding  chapters. 

In  Chapter  III  the  mathematical  theory  concerning  the  solution  for 
the  optimum  driving  function  for  the  linear  portion  of  the  relay 
servo  loop  will  be’ treated  in  detail.  Chapter  IV  will  then  describe 
the  digital  computer  programs  used.  Chapter  V  will  present  the  quan¬ 
titative  analysis  of  the  Honeywell  controller  by  the  methods:  of 
Chapter  III*  An  analogue  computer  simulation  and  its  results  will 
be  set  forth  in  Chapter  VI  to  form  a  basis  for  comparison  of  the 
results  of  the  preceding  section.  Comparison  of  these  pesults  and 
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conclusions  will  be  left  to  Chapter  VII.  Related  to  . Chapters  III  and 
IV  will  be  three  appendices.  Appendix  A  will  explain  in  detail  the 
various  FORTRAN  statements  of  the  digital  computer  programs;  in 
Appendix  B  comparative  solutions  of  an  example  problem  simple  enough 
to  be  worked  out  by  hand  will  be  included;  Appendix  C  will  show  the 
solution  of  a  typical  fourth  order  problem.  The  numerical  results 
of  these  last  two  appendices  will  demonstrate  to  the  reader  t'le  utility 
and  validity  of  the  methods  presented.  Tabular  data  will  be  presented 
in  Appendix  D. 
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II,  Description  of  the  Minneapolis-Honevwell  Adaptive  System 

Before  analysis  of  the  Minneapolis-Honeywell  Adaptive  Control 
System  by  State  Space  techniques  it  would  be  well  to  consider  the 
system  in  a  qualitative  manner  so  that  the  application  of  State  Space 
analysis  will  be  clear. 

It  should  be  noted  here  that  important  symbols  used  throughout 
this  report  are  listed  and  defined  in  the  Table  of  Symbols  in  the 
prefatory  part.  Figure-  3  on  the  following  page  is  a  block  diagram 
of  the  Honeywell  pitch  rate  system  which  was  installed  in  the  F94C, 
This  block  diagram  was  taken  directly  from  the  WADC  report  (Ref  B:85), 
The  transfer  functions  of  the  various  components  are  shown  in  their 
respective  block's.  These  correspond  to  Flight  Condition  13  in  the 
WADC  report  (Ref  8:153).  Flight  Condition  13  refers  to  the  F94C 
aircraft  operating  at  an  altitude  of  35,000  feet  and  at  a  mack  number 
of  .73.  Note  that  ©  as  used  in  this  chapter  refers  to  pitch  rate. 

The  operation  of  the  pitch  rate  system  is  discussed  in  the  following 
section. . 

Operation  of  the  Adaptive  System 

*he  basis  for  operation  of  the  adaptive  control  system  shown  in 
Figure  3  in  the  forcing  of  the  aircraft  pitch  rate  response  Q  to 
follow  the  response  of  the  model  as  closely  as  possible.  The 
model  is  a  simple  electrical  network  with  response  characteristics 
equal  to  those  a  pilot  would  desire  over  the  entire  flight  envelope 
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of  the  aircraft,  A  second  order  function  with  a  damping  ratio  of 
.63  and  a  natural  frequency  of  3.17  is  used.  The  desired  model  response 
is  fed  into  the  relay  servo  loop.  The  term  relay  servo  loop  here 
refers  to  that  portion  of  the  system  enclosed  by  the  rate  gyro  feedback, 
and  includes  tha  rate  gyro  itself.  The  model  output  is  summed  with 
the  negative  .of  the  aircraft  pitch  rate  which  is  fed  back  by  means 
of  the  rate  gyro.  The  rate  gyro’s  characteristics  are  such  that  this 
is  essentially  unity  feedback.  An  error  signal  is  thus  produced. 

The  error  signal  is  fed  through  the  switching  logic,  which,  as  will 
be  shown  in  Chapter  V,  is  approximately  a  proportional  plus  derivative 
network.  The  algebraic  sum  of  the  error  signal  and  of  its  derivative 
multiplied  by  ,225  controls,  the  switching  of  a  very  sensitive  electronic 
relay.  If  the  sum  is  positive  the  output  of  the  relay  will  be  positive, 

A  2000  cycle/second  sinusoidal  dither  signal  is  applied  to  the  relay 
in  order  to  improve  its  characteristics.  The  gain  changer  and  the 
limiter  provide  a  means  for  cutting  down  the  relay  output  voltage 
for  small  error  signals.  In  addition  thp  limiter  assures  proper 
limiting  of  the  output  of  the  electronic  relay. 

The  filter,  the  servo  and  actuator,  and  the  aircraft  compose 
the  linear  controller  portion  of  the  system.  Transfer  functions  shown 
for  the  servo  and  actuator  and  for  the  aircraft  were  derived  experimentally. 
The  filter  is  a  lag-lead  network  used  to  improve  the  response  of 
the  linear  system.  The  denominator  provides  a  pole  close  to  the 
zero  of  the  servo  and  actuator  thus  preventing  this  zero  from  greatly 
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affecting  the  response  of  the  loopf  the  numerator  assures  that 
roots  due  to  the  servo  and  actuator  poles  do  not  become  dominant. 

The  relay  servo  loop  of  the  Honeywell  adaptive  system  closely 
resembles  that  of  the  practical  adaptive  system  depicted  in  Figure  2 
in  the  introductory  chapter.  The  differences  are  the  added  components 
of  the  Honeywell  loop;  these  are  the  limiter,  the  gain  changer,  and 
the  AC  dither.  These  components  provide  second  order  improvements  to 
the  system  and  their  inclusion  complicates  the  analysis  unnecessarily. 
Thus  they  will- not  be  included  in  any  further  analysis  of  the  system. 

Also,  the  response  of  the  model  can  easily  be  obtained  by  the 
Laplace  transform  technique,  or  any  one  of  several  others.  Therefore, 
only  the  relay  servo  loop  will  be  analyzed  in  this  report. 

Before  going  on  to  Chapter  III  the  reader  should  note  that  the 
filter,  the  actuator  and  servo,  and  the  aircraft  comprise  the  linear 
controller  for  which  an  optimum  driving  function  is  desired.  The 
next  chapter  of  this  report  will  develop  the  mathematical  theory  for 
finding  this  optimum  driving  function. 
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III.  Mathematical  Development  of  the  Solution  for  the  Optimum  Driving  Function 

Since  the  use  of  State  Space  techniques  and  the  solution  for 
the  optimum  driving  function  are  both  rcJ.atively  new  this  portion 
of  the  report  will  be  devoted  to  a  detailed  development  of  the  under¬ 
lying  mathematics  covering  this  problem.  Because  matrix  notation  will 
be  used  freely  throughout  the  development  it  might  be  well  to  define 
the  symbols  to  be  used.  A  matrix  will  be  represented  by  an  underlined 
upper  case  letter,  i.  e.,  ^  .A  vector  (n  x  1  or  1  x  n  matrix)  will 
be.  denoted  by  an  underlined  small  letter,  i.  e.,  .  A  transposed 

matrix  will  be  denoted  by  ^ ,  as  0^.  The  symbol  ^  will  be  used  to 
signify  the  inverse  of  a  matrix,  i.  e.,  0  1  Successive  time  derivatives 
will  be  denoted  by  successive  dots  above  the  symbol  such  as  A.  or  (J  . 

The  familiar  derivative  notation  ^  and  the  operator  D  will  also 
be  used.  The  Laplace  operator  will  be  denoted  by  5  •  Vector  or  matrix 
elements  will  often  be  written  in  subscripted  form  for  clarity  of 
explanation.  Lower  case  letters  will  always  be  used  in  this  instance. 

The  first  step  toward  finding  the  optimum  driving  function  is 
the  development  of  a  transformation  which  will  make  it  possible 
to  put  a  general  linear  transfer  function  into  State  Space  notation. 

This  transformation  will  be  developed  in  the  next  section  of  this 
chapter. 

Transformation  of  a  General  Linear  Transfer  Function  into  State  Space  Form 

In  the  introductory  chapter  ?  transformation  was'  introduced  to 
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allow  an  Nth  order  differential  equation  to  be  written  as  a  series 
of  linear  first  order  differential  equations.  However,  it  often  occurs 
in  control  system  problems  that  the  right  hand  side  of  the  differential 
equation  contains  derivatives  of  the  driving  function.  For  instance 
the  transfer  function 
gp  _  5  +  eL. 

becomes 

+  ce.  a-e^  . 

Since  derivatives  of  many  inputs,  for  instance  step  or'  impulse  functions, 
are  not  calculable  in  the  time  domain  it  is  necessary  to  replace 
the  (differential  system  by  a  first  order  system  containing  no  derivatives 
on  the  right  hand  side. 

A  transformation  has  been  developed  by  Lanning  and  Battin  (Ref  6:191} 
for  a  general  time  varying  system  to  eliminate  right  hand  side  derivatives. 
Since  this  paper  deals  only  with  the  time  invariant  or  constant 
coefficient  case,  a  formula  will  be  developed  inductively  considering 
the  coefficients  of  the  differential-  equation  to  be  constant.  Consider 
the  differential  equation 

The  following  system  of  first  order  equations  is  then  to  be  made 
equivalent  to  equation  (16): 

7-  -  "X,  +  Qro'U'  (17) 
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I 


t  Cr, 


Eliminate  X,  and 


-  ^  j  I  +  6'i, 

successively  in  the  following  manners 


(18) 


(19) 


7(;  -  X-  GroXC 

Differentiate  (20)  with  respect  to  times 


r  ^ 


Substituting  this  result  in  (18)  and  rearranging 


(20) 


(21) 


■Jf  Or,  AA^ 


(22) 


Now  differentiate  (22) s 


After  rearranging  (23)  substitute  the  result  for 
ranging  the  resulting  expression  then  yields 


(23) 


n  (19)»  Rear- 


~  +  Crp'Cc,  (24) 
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f 


r 


Now  eliminating  Xf  and  Xi  by  substituting  (20)  and  (22)  one  has 


t  (I, 


+  eCi,x  r  Cr, 


^  G^.A^  (25) 


■>c  +  Cc, 


Equating  right  hand  coefficients  of  (25)  with  those  of  (16)  yields 


ba  =  &-0 

..  h,  ~  G■^  t  Cro 
bi  =  6-2,  +-  G-o 

Solving  for  (7, ,  Gr,  *  and  (y^  then  yields  ■ 

60  = 

G-|  -  -  A,  6-0 

Gr  j  ~  ^  1 6-0 

Inductively,  the  following  relation  can  then  be  derived: 
i-  I 

'  C--  ;  ...  .n 

k-o  ^  > 


(26) 

(27) 

(28) 

(29) 

(30) 

(31) 

(32) 


Thus  the  N  order  system 


(33) 
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may  be  transformed  into  a  series  of  first  order  differential  equations 

of  the  form 

(34) 

1 

X,  -  TC.t.'h  0-^  yCc- 

4 

(35) 

t 

(36) 

where  &.C  (-t*  0,  lj> • 'j/W )  is  given  by  formula  (32)  and /H 

is  the  order 

of  the  system.  In  matrix  form  this  is  written  as 

i.(^)  =  -h 

(37) 

j^io)  =  C 

(38) 

where  the  state  vector  and  ■)c,,:c^««»^^are  the  state  variables. 

So  far  a  method  has  been  presented  for  expressing  a  general 
linear  transfer  function  in  State  Space  notation.  It  is  now  necessary 
to  consider  this  transformed  transfer  function  in  its  place  as  the 
linear  controller  of  a  relay  servo  loop.  As  was  pointed  out  in  the 
introduction,  it  is  presently  impossible  to  synthesize  a  circuit 
which  would  perform  the  switching  logic  function  in  an  optimum  manner. 
What  can  be  done,  however,  is  the  following;  given  a  linear  controller 
as  described  mathematically  by  equations  (37)  and  (38)  it  is  possible  to. 
compute  the,  driving  f unction/u^which  will  reduce  the  state  vector  'X(t) 
to  zero  in  the  shortest  possible  time.  If  one  'lets  the  state  variables 
represent  the  error  of  the  system'^  then  one  will  have  computed  the 
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driving  function  for  the  linear  controller  which  will  reduce  the 
error  to  zero  in  the  shortest  possible  time. 

Thus,  although  an  optimum  switching  logic  cannot  be  synthesized 
the  output  of  such  a  switching  logic  can  be  found.  Finding  this 
optimum  driving  function  is  a  useful  tool  of  analysis  since  the  output 
of  the  controller  when  subjected  to  this  driving  function  can  be  found. 
Comparisons  can  then  be  made  between  these  inputs  and  outputs  and 
those  of  the  same  linear  controller  when  utilized  in  a  practical 
system. 

Calculation  of  the  Optimum  Driving  Function 

The  problem  to  be  considered  in  this  section  is  the  calculation 
of  the  optimum  driving  function  for  the  system  of  .equations  (37) 
and  (38).  A  method  developed  by  Ho  and  Brentani  was  followed  to  solve 
this  problem  (Ref  5)*  First,  consider  a  mathematical  statement  of  the 
problem. 

Given  a  linear  system 

=  f-KU)  -f  2 

-x[0)  -  Q  (-1°) 

where  C  is  an  initial  error  in  the  system,  compute  the  driving  function 
necessary  to  minimize  the  error  in  an  optimum  manner.  It  will 

be  necessary  to  state  mathematically  the  criterion  or  norm  used  in 
minimizing  the  error.  However,  this  statement  of  the  criterion  function 
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will  be  made  in  a  more  appropriate  part  of  this  discussion<  The  general 
method  of  solution  for  the  problem  must  be  discussed  first* 

The '  l^umerical  Approximation*  Since  a  closed  solution  to  this 
problem  is  in  the  general  case  presently  unavailable,  a  numerical 
approximation  will  be  made  and  a  digital  computer  will  be  employed 
in  the  computation.  Therefore,  it  will  be  necessary  to  divide  the 
control  interval  into  a  number  of  subintervals  and  make  the  driving 
function piecewise  constant.  Thus 


X* 


T-0 

K 


A  r  A  ;fc  <  ;  A  -  K~  I 

(42) 


where  "C  is  the  step  length. and iK  the  number  of  control  ^subinterwals.  One  must 
specify  the  control  interval  T  to  be  used  for  each  solution,  the 
control  interval  may  also  be  referred  to  as  the  terminal  time  in 
this  report. 

Each  now  becomes  a  variable  in  itself  and  the  problem  becomes 
one  of  many  variables,  and  of  minimizing  a  function  of  these  variables. 

To  do  this,  a  method  of  successive  approximations  known  as  the  Steepest 
Descent  method  is  employed. 

The  Steepest  Descent  Method.  Before  entering  into  a  mathematical 
discussion  of  the  method  of  Steepest  Descent  an  intuitive  approach 
should  prove  helpful.  Figure  4  on  the  next  page  is  a  map  of  contour 
lines  of  a  three  dimensional  surface'  which  is  shaped  like  a  bowl. 

If  one  were  at  a  point  on  the  rim  and  wished  to  get  to  the  lowest 
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point  as  fast' as  possible  it  would  be  better  to  go  straight  down  the 
slope  in  the  direction  of  the  center  arrow  than  to  spiral  around 
the  rim  in  the  direction  of  one  of  the  other  two  arrows.  This  is 
assuming  that  there  aren't  too  many  obstacles  blocking  the  path  of 
descent.  It  may  be  deduced  even  further  that  even  if  one  cannot  see 
the  bottom  of  the  bowl  one  could  continually  pick  the  steepest  slope 
and  go  down  this  slope  until  he  gets  to  the  bottom,  that  is  until  he 
starts  going  up  again.  The  direction  taken  in  this  example  is  the 
negative  of  the  gradient. 

The  problem  may  now  be  restated  mathematically  in  a  more  general 
form.  Consider  a  function  to  be  describable  in  dimensional  Euclidian 
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space.  One  then  computes  the  gradient  of  the  function  at  some  point 
on  a  hypersurface  in  that  space.  Descending  from  that  point  in  space 
along  the  negative  of  the  gradient  then  gives  the  greatest  rate  of 
decrease  of  the  function.  Although  this  follows  from  the  previous 
discussion,  a  mathematical  proof  is  in  order (Ref  7:2-4)»  Consider 
a  function  df  several  variables  ...  Start  at  some  point  Xtio) 
in  /*\  dimensional  space  and  move  an  infinitesmal  distance  t 

The  total  differential  of  Xii)is 

JJ  -  AJl 

^  >1  3  ^  “Xwf 

From  the  above  the  derivative  with  respect  to  ^  can  be  expressed  as 


/A 


(45) 


^  =  9  U  . 

One  desires  to  find  the  most  negative  rate  of  descent  ^  subject  to 


the  constraint  of  (43).  In  order  to  do  this  first  rewrite  (43)  as 


(46) 


Then  adjoin  (46)  to  (45)  by  means  of  a  Lagrangian  multiplier (Ref  4;120-125); 


iH.’k  - 


(47) 


»• 


r 


V 
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Take  partial  derivatives  in  (47)  with  respect  to 


^^and  set  each 


equal  to  zero: 


The  following  is  obtained  by  rearranging  (40): 


(49) 


Now  form  the  sum  of  the  squares  of  each  of  the  equations  of  (49): 


(50) 


The  left  hand  side  of  (50)  is  seen  to  be  unity  from  equation  (43): 


(51) 


Solving  (51)  for  yields 


This  value  for  X  can  be  substituted  into  (49)  to  obtain 


(53) 


Recalling  that 
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Substituting  (45)  into  (53)  after  solving  for yields  the  most 
negative 


(54) 


The  inner  summation  is  performed  before  the  outer  summation  and  is 
therefore  factorable: 


Again  recalling  that 


The  second .summation  becomes  /  and  thus 


(55) 


(43) 


Examination  of  (56)  reveals  that  these  are  respectively  the  positive 
and  negative  gradients.  The  maximum  rate  of  descent,  which  is  the 
maximum  rate  of  change  of  the  function  with  respect  to  distance,  is 
therefore  obtained  by  making  the  direction  of  descent  equal  to  the 
negative  of  the  gradient. 


Complete  Statement  of  the  Control  System  Problem*  With  the  above 
result  in  mind  it  is  now  possible  to  consider  the  problem  of  obtaining 
a  solution  for  the  optimum  driving  function  of  the  linear  controller 
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by  the  Steepest  Descent  method*  In  an  actual  autopilot  system,  an 
additional  constraint  must  be  added  to  those  already  stipulated  in 
the  preceeding  sections.  This  additional  constraint  limits  the  maximum 
excursion  of  the  driving  function  to  one  that  will  not  overdrive  the 
system,  A  more  complete  mathematical  statement  of  the  problem  may 
now  be  made.  Consider  the  dynamic  system 


(57) 

Z(P)  '  1 

(50) 

which  can  be  written  in  subscripted  form  as 

/V 


(59) 

(60) 


It  is  now  necessary  to  state  the  function  to  be  minimired: 

V=  II  II  ^  ^  -  ’XlT)  B-  79(J)  (61) 

wherd  ^  id  a.  positive. definite  symnetric^l  matrix  of  weighting  functions. 
The  functioh  (61.) , will  be  referrsd.tb  as  the  parformance  function 
or  as  the  terminal  norm.  Since  V  is  quadratic  in  form,  finding  the 
minimum  V  will  also  find  the  minimum  of  X(T)»  The  choice  of  the  matrix  S, 
will  be  discussed  in  detail  in  Chapter  IV, 

The  system  will  be  subject  to  the  maximum  input  constraint 
mentioned  above: 

Iac(^]I  Y  (62) 
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or  for  Afy  piecewise  constant 


6  Y 


As  Oj  tj 


) 


K-l 


(63) 


Determining  the  Gradient*  Now  that  a  complete  mathematical  statement 
of  the  problem  has  been  made  the  method  of  solution  may. be  considered. 
Briefly  stated,  the  method  of  solution  is  to  compute  the  gradient 
of  V  wiih  respect  to  and  then  to  descend  along  this  path  subject 

to  the  constraint  that  may  not  go  beyond  a  boundary  foamed  by 

equation  (63). 

As  before  define 


ir  Iro 


4tz  ^ x.<  \)  t 

K-i 


(41) 

(42) 


Thus 


- 


Ay, 


(64) 


It  is  necessary  to  find  the  change  in  the  terminal  norm  V  with  respect 
to  a  change  in  Recalling  that  each>Mi|is  a  separate  variable 

in  K  dimerfsional  space  the  following  relationship  can  be  formed: 


iV  .  J’lir'  hi(T;  A  o-;K-l 

.This  says  that  for  any the  change  in  with  respect  to  that 


(65) 
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equals  the  change  in  V  with  rsapact  to  the  change  in  the  terminal 
state  vector  multiplied  by  the  change  in  the  terminal  state  vector 
with  respect  to  the  change  in.  the  control  •  At  this  point  define 


^^iCT)  _  ^ 

as  an  element  of  an x  K  matrix  H«  Making  this  substitution  in 
(65)  yields 


A-jij 


Oy,K-( 


(67) 


Performing  the  indicated  summation 


iC  ' 


In  matrix  form  this  is  written  as 


k-  K-l 


dV  .  M^av^ 
a^.  .  '  ^  hZ:CT) 

Recalling  the  expression  for  the  tenninal  norm 

V  =  S  ZS- 


(69) 


(61) 


For  R  symmetrical-  the  derivative  is  calculated  to  be  (Ref  3:4B) 

^  a.  R  2C*; 

— 

Substituting  (70)  into  (69)  yields 


(70) 


4^  ^  ^h'i^xLt) 


(71) 
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X(T^can  be  calculated  by  solving  the  system  aquations*  Thus  it  is 
only  necessary  to  determine  ^  in  order  to  find 


Calculation  of  fcl  .  An  examination  of  the  eolution  to:, the  system 
equations  will  lead  to. e  method  for  calculating  h!  •  solution  to 
equations  (57)  and  (50)  is  (Ref  lt78) 

2(r)  =  I  (Tjo;  c  +  Jj  I  cz  jg-  (72) 

where 

i  - 1  (73) 

where  X  is  the  identity  matrix*  To  find  $^l^which  is  under  the  integral 
sign  the  adjoint  system  may  be  utilized*  This  is  defined  as 

-£i(r,<r)  >  |«--x 

It  is  now  necessary  to  consider. the  solutions  to  equation  (73)  and 
to  the  adjoint  aquation  (74)*  Equation  (73)  resembles  the  scalar  form. 


The  solution  of  the  above  scalar  equation  is  known  from  differential 
equation  theory  to  be 


at 


where  the  exponential  is  defined 


(76) 
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In  like  -manner  the  solution  to  equation  (73)  is 

^Ct)  ~  X 

if  is  defined  as  follows 

1  ;  £"  X  C79> 

yti:c>  yh{  - 

That  this  series  converges  and  is  differentiable  term  by  term  is 
readily  shown  (Ref  ,3:42-45)»  Term  by  term  differentiation  and  sub¬ 
stitution  into  equation  (73)  shows  that  equation  (79)  is  indeed  a 
solution!  The'  matrix  $  is  "often  known  as  a  fundamental  matrix  for 
the  system.  Similarly  the  solution  to  (74)  is  found  to  be 

§CT-(r)  r  g_ECT-tr)  (eo) 

Substitution  of  these  results  into  (72)  yields 

Equation  (01)  is  the  same  as  the  solution  (72)  except  that  it  is 
written  in  exponential  foim»  The  definitions  of  the  exponentials 
may  be  used  to  calculate  a  numerical  solution  for  7^*  However,  this 
problem  will  be  taken  up  in  Chapter  IV,  Use  of  the  exponential  form 
will  also  clarify  the  following  steps  in  the  development  of  a  method 
for  the  calculation  of  the  matrix  ^  » 

Suppose  for  a  moment  that/<c{<^wsre  constant  over  the  entire  control 
interval.  It  then  could  be  factored  from  under  the  integral  thus  facilitating 
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differentiation  with  respect  to  AL,  • 

Z(-T)  =  ^  ^  oC(r 

The  derivative  for  >tc  constant' now  becomes  ( 


I 


(03) 


The  value  of  making  Aa^  piecewise  constant  now  becomes  evident. 

The  integration  can  he  broken  up,.:  wi-jih  the  intetval  of  each  part 
chosen  to  be  the  control  subinterval  T  p  and  each  ^4_  can  be  factored 
from  the  integral.  If  the  exponentials  are  chosen  judiciously  the 
following  result  is  obtained: 


XCt)  = 


-f 


t 


X<r  (64) 

■  ’ll ^  tiff 


Now  differentiate  XC<1  with  respect  to  each  : 


i_CjCr)l.  ^  cc«.w;r-(r3/  ^  (sa) 

The  exponential  is  a  constant  and  may  be  brought  from 

under  the  integral: 


f  -tr]f  ^  ^  (b6) 


The  integral  of  (66)  can  be  simplifisd  to  a  much  more  maneagable 
form.  To  do  this  it  is  first  necessary  to  evaluate  the  integral* 
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To  do  this  perform  the  integration; 

-T)  1 

By  equating  this  result  with  that  obtained  from  evaluation  of  the 
following  integral  a  convenient  way  of  calculating  the  integral  of 
equation  (87)  can  be  determined; 


(87) 


Integration  yields 

i  =  »’> 

Equating  the  results  of  the  two  integrations  one  obtains 

^  Ur  .  f 

Substituting  the  results  of  (90)  back  into  equation  (86)  then  yields 

— 

Reverting  back  to  the  notation  of  the  fundamental  matrix  one  has  the 
following  which  is  equivalent  to  (91): 

}  (r~  ^  f.  JiT 

d  ^4.  "■  )o  ^ 


(92) 
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In  order  to  simplify  the  notation  let 

2:= 

/fl 

Recalling  that 

Putting  (92)  into  matrix  notation  one  has  the  result 

^  ~  ^  (T~^Z)  A  .  »  *  (94) 

Distance  to  Descenda  Having  determined  a  means  for  calculating 
the  gradient  it  is  now  necassary  to  determine  the  proper  distance 
to  descend  along  this  gradient  in. order  to  maximize  the  change  in 
the  terminal  norm  V  •  Putting  the  solution  to  the  system  equations 
in  a  more  maneagable  form  facilitates  the  algebra  involved  in  finding 
this  proper  distance  to  descend.  Recall  the  solution  to  the  system 
equations: 

'^Lr)  -  ey  ^{<r)  OD 

Let 

P  ^  e^'^C  (95) 

Then  using  the  result  for  f~^  ,  which  is  equation  (94)  and  the  matrix 
version  of  the  driving  function  which  is  equation  (64),  a  matrix 
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formula  for  'iC(j)caT\  be  obtained: 


%(.-r)  -  P  + 


(96) 


Now  to  continue  with  the  calculation  of  the  proper  distance  to  descend 
along  the  gradient  consider  the  terminal  norm  formula,  which  is  equation 
(6l),  with  equation  (.95)  substituted: 


'/=  (£'  f  R  ) 


(97) 


However,  the  change  in  the  terminal  norm  v  with  a  change  in  the 
driving  function  is  of  interest.  Taking  ✓<<.  at  any  point,  say  ■C' , 
and  performing  the  matrix  operations  indicated  in  (97) 


=  p'R  p  +  (98) 


Consider  a  small  variation  in 

+■  (99) 


Substituting  (99)  into  (98)  yields 

Vet,’  f'Re  V  *p'RH(^ci4^iDo) 

Then  carry  out  the  multiplication  called  for  in  (100) 

vi.+/  =  +  ^2  M'B  8  ■f'  z  p'  R  8 

H' R  8  .+  x^'RpSiy 
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Now  define 

SV-  V<,^I~  (102) 

Then  subtracting  (90)  from  (101)  leaves 

SI/--  4  H  S±r  (103) 

Combining  ! the  first  twb  teiws 

Sl/=  -I-  H' 

Substituting  equation  (96)  into  (104) 

I 

S^- )£H  5i«-  +  $iJB'SdS-s- 

Recalling  equation  (71)  and  taking  its  transpose  one  has 

Substituting  the  above  result  into  (105) 

^1/  --  l^ySA^-  + 

Now  define 

^  (106) 

Finally 

Si'  -  +  Sjf-' H' a  B  <“’> 
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Equation  (109)  is  an  expression  for  the  variation  in  the  terminal 
norm  vetses  a  small  variation  in  the  driving  function  It  is  now 

desirable  to  maximize  the  change  in  the  terminal  normal' by  determining 
a  proper  distance  to  descend  along  the  gradient  The  change  in 

the  driving  function  equal  to  the  direction  of  descent  multiplied  by 

the  proper  distance  to  descend) 


-  'y}  ^  (iio) 

where  yf  is;  a  constant  and  is  the  proper  distance  to  descend* 

Were  it  passible  to  descend  along  the  gradient  without  being 
subject  to  the  constraint  on  the  driving  function  discussed  on  page 
25,  yj  could  be  determined  in  the  following  manner.  Substitute  (110) 
into  (109): 


sv^ 

Calculus  can  now  be  employed  to  determine  the  value  for  7)  which  will 
maximize  gl/*  Differentiate  equation  (111)  with  respect  to  >7  and  set 
the  resulting  expression  equal  to  zero< 


r  f'(h  -y  ^  ti'  KH  ^  ^  ^ 

Solving  for 

'Vf  ^ 


(112) 


(113) 
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However,  limited  since  ^y£^ust  not  go  outside  the  boundary 

which  was  determined  by  the  maximum  contraint  on  the  driving  function  <ci>. 
Proper  use  of  this  constraint  will  aid  in  computing  the  proper  change 
in  the  driving  function  subject  to  the  constraint* 

The  Descent  Sehynas*  Two  different  descent  schemes  may  be  de¬ 
termined  utiliring  i|t,he  constraint  on  the  driving  function.  These  are 
the  Gradient  Projection  Descent  Scheme  and  the  Corner  Aiming  Descent 
Scheme.  Both  schemes  use  the  constraint  on  the  driving  function  to 
limit  hut  they  do  it  in  a  slightly  different  way.  A  geometrical 

I 

interpretation  of  the  Gradient  Projection  Descent  Scheme  is  shown  in 
Figure  4a  which  is  is  two  dimensional  model  used  to  aid  in  the  ex¬ 
planation  of  the  descent  schemes.  Figure  4a  is  located  on  the  next 
page.  Since  this  is  a  two  dimensional  modal  only  two  control  sub- 

'  I  '  ■ 

intervals  are  used.!  As  before,  the  attempt  iS' made  to  mihimize  the 

i  < 

terminal  norm  V  byj  changing  the  driving  function  The  negative  of 

the  gradient  is  thelrefore  drawn  from  point  A  across  the  lines  of 

■j  ■  ' 

constant  \/  which  airs  the  same  as  the  contour  lines  of  Figure  4. 

•I  '  '  i 

I  .  ' 

The  vector  AB  thus  jcepresents  <^with  7^  computed  as  in  squstion  (113). 
However,  the  boundary  For  the  driving  function  does  not  allow  a 

I  .  ' 

f  •  *  I  j 

vector  of  this  lehg|th.  The, tip  of  the  vector  is  therefore  projected 

down  to  the  boundjsr^  to  point  C.  A  vector  drawn  from  point  A  to 

! 

the  intersection  of .  the  projection,  C,  is -’then  used  as  the  change  in 
control.  It  is  frbmi  this  procedure  that  the  term  "gradient  projection" 
is  derived. 
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Notice  that  a  new  gradient  from  point  C  can  be  determined  and 
the  tip  of  the  vector  CO  projected  to  the  boundary  as  before.  This 
time  movement  from  C  to  E  is  along  the  boundary.' This  movement  along 
the  boundary  leads  one  to  the  conclusion  that  for  the  orientation 
of  the  control  boundaries  oh  the  contour  map  of  Figure  4a,  the  driving 
function  will  always  go  to  the  corner  of  the  boundaries.  The  termination 
of  the  driving  function  44  at  the  corner  of  the  boundaries  is  in 
keeping  with  the  maximum  principal  of  Pontriagin  which  states  tha-t 
the  magnitude  of  the  driving  function  ^^must  be  the  maximum  allowed 
by  the  physical  limitations  of  the  system  in  order  to  achieve  time 
optimization.  A  second  deecent  scheme,  the  Corner  Aiming  Descent  Scheme, 
can  be  formulated  by  taking  advantage  of  this  principal.  Figure  4b 
illustrates  this  second  descent  scheme  geometrically.  ' 


Geometrical  Interpretation  of  Comer  Aiming  Descent 
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After  the  gradient  is  determined,  instead  of  projecting  the 
tip  of  the  vector  down  to  the  boundary,  the  corner  nearest  the 
vector  AB  is  aimed  for.  In  this  way,  the  assumption  that  the  corner 
will  eventually  be  reached  is  utilized  to  aid  in  the  descent.  The  use 
of  the-  maximum  principal  thus  enters  into  the  calculations. 

All  of  the  statements  mads  about  the  two  dimensional  models  of 
Figures  4a  and  4b  can  of  course  be  extended  to  K  dimensional  hyperspace 
for  the  K  number  of  control  subintervals  encountered  in  an  actual 
problem,  ^he  mathematics  for  computing  the  change  in  the  driving 
function  utilizing  these  two  descent  schemes  is  presented  below. 

The  following  equations  describe  the  Gradient  Projection  Descent 
Scheme: 


Yjt'/CCn  <  (114) 

-I 

where  Vyf  'the  upper  boundary  for  the  driving  function  and  ^  is 
the  lower  boundary  for  the  driving  function.  A  new  yj  must  now  be 
determined  in  the  same  manner  as  before: 

V'  -  — 

It  is  now  necessary  to  assure  that  does  not  exceed  unity; 
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Therefor* 


where 

7-  /^/^/ 

r  i-/  IF  IX)>  I 

Finally,  the  change  in  the  driving  function  is  expressed  as 


(116) 


(117) 


Conceptually,  the  above  procedure  limits  as  shown  in  Figure  5,  below* 
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The  following  equations  formulate  the  Corner  Aiming  Descent  Scheme: 


'«rt  < 


"=  ° 


As  before  it  is  necessary  to  determine  a  new 


B' Bh  ^ 


Again  must  be  limited  to  unityl 


(119) 


(120) 


(121) 


Finally t  as  before 

(122) 

In,  the  Corner  Aiming. Descent  Scheme  u/" is' limited  in  a.  different' 
i^ay.  Figure  6  on  .  the  next  page;  shpwf,  this.  concaptua4.1yy;  : 

The  new  value  of  thf;  driving  functipn  csn.  now  be  computed  by 
adding  the  value  for  5^  computed  by  either  the  Gradient  Projection 
Descent  Schema  or  thahCpmer'-Aiming  , Descent  Scheme 
to  the  old  value  of  ^ •  This  is  expressed  mathematically 

~  S (123) 
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♦ 


Figure  6 

Corner  Aiming  Descent 

Completion  of  Solution 

The  process  is  now  repeated  with  the  new  value  of  and  the 
solution  will  hopefully  converge  to  a  true  minimum.  It  is  necessary^ 
however,  to  determine  when  the  minimum  has  been  reached.  In  the 
method  of  Steepest  Descent  a  minimum  is  reached  as  soon  as  the  gradient 
becomes  zero  as  in  Figure  6a  on  the  next  page.  However,  in  this 
problem  the  constraint  on  the  driving  function  may  prevent  all  the 
elements  of  the  gradient  vector  from  reaching ^tero.  Therefore,  conditions 
other  than  a  zero  gradient  will  have  to  be  used  to  determine  when 
the  minimum  for  the  terminal  norm  V  has  been  reached.  The  solution 
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is  considered  complete  whenever  one  of  the  three  conditions  below 
occurs  for  all  /U*: 


-  0 

If 

<0 

. 

Kt  <  yUA 

If 

0 

II 

(124) 

-yCi^  =  0 

IF 

^4  >  0 

Equation  (124)  states  that  the  solution  is  complete  when,  for 
each  control  subinterval,  l)  the  gradient  has  reached  zero  and  no 
further  improvement  in  the  terminal  norm  V  can  be  made,  or  2)  that 
no  further  change  is  possible  in  the  driving  function  .<<4, due  to  the 


constraint  imposed.  When;. one  of  these  conditions '  is  met  for  every 
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control  subinterval  computation  is  halted  and  the  final  driving 
function  is  read  out*  That  the  point  reached  is  a  minimum  is  shown 
in  Reference  5  (Ref  S:A-2DtoA^23)» 

It  should  be  noted  that  any  time  the  method  of  Steepest  Descent 
is  employed  the  minimum  reached  may  only  be  a  local  minimum.  There 
may  still  exist  someplace  on  the  hypersurface  considered  a  global 
minimum.  Figure  6a  illustrates  this  for  a  two  dimensional  problem. 

This  should  not  detract  from  the  usefulness  of  this  method,  however, 
since,  for  the  type  of  problem  being  considered  in  this  report,  the 
approximate  -solution  is  already  known  from  the  maximum  principal* 
According  to  the  maximum ' principal  the  solution  is  expected  to  be 
"bang-bang"  and  thus  an  approximate  check  on  the  results  is  available* 
This  completes  the  mathematical  development  for  the  solution 
for  the  optimum  driving  function  of  a  linear  controller*  The  next 
chapter  of  this  report  will  present  a  description  of  the  computer 
programs  used  and  a  general  discussion  of  their  use. 
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IV.  Description  of  the  Computer  Programs 


The- purpose  of  this  chapter  is  to  extract  from  Chapter  HI 
the  various  equations  utilized  in  the  digital  computer  programs 
and  present  them  in  a  meaningful  manner j  in  addition,  certain  here¬ 
tofore  unexplained  phhses  of  the  solution  will  be  treated  in  detail. 
Since  the  development  of  the  equations  of  the  descent  has  been  set 
forth  in  great  detail  in  the  preceding  chapter,  the  equations  extracted 
from  that  chapter  will  be  stated  without  further  justification. 

Appendix  A  will,  by  reference  to  the  listings  of  the  programs,  explain 
the  purpose  of  the  various  FORTRAN  statements  used. 

Computation  of.^  «  flr  .  and  H 

Although  p  ,  and  H  were  defined  in  Chapter  III  the 

details  of  their  computation  remain  to  be  discussed.  Recalling  the 
definitions  of  p  ,  g,  ,  and  H 


(95) 


«■  --  I  (.r)  J  <”> 

-  "  (T-r)a  ^Cr-i-c)  ^  (94) 


The  most  obvious  way  to  compute  these  functions  is  to  use  the  definition 


of  I  ; 

f  u;  =  \  =  I 


(78) 
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Equation  (TB)  can  easily  be  calculated  by  usings  the  definition  of 


(79) 


The  successive  terms  of  the  series  are  computed  and  summed  on  a 
digital  computer  for  a  particular  value  of  This  method  was 
originally  used  for  computing  ^  and  worked  quite  well  for  problems 
with  a  short  control  interval  and  small  moduli  of  the  elements  in 
the  matrix  F  .  Due  to  the  arithmetical  nature  of  the  digital  computer 
it  is  necessary  to  solve  the  numerical  rather  than  the  algebraic 
problem.  For  instance,  for  each  column  of  tl  a  separate  ^fr-4t'3must 
be  computed* 

The  number  of  terms  in  the  series  for  ^  necessafy  for  the 
required  accuracy  varies  with  the  moduli  of  the  elements  in  the 
numerator  of  the  terms  of  the  series.  (Note  that  for  the  actual 
computation  the  matrix  P  is  multiplied  through  by  the  scalar 
before  the  numerator  is  carried  to  a  power. )’ The  size  of  the  . 
numerator  moduli  depends  not  only  on  the  elements  of  the  matrix 
but  also  on  the  value  chosen  for  Experience  has  shown  that  whenever 
it  is  necessary  to  carry  the  series  farther  than  100  terms  the 
computation  time  is  .excessive  and  roundoff  errors  in  the  computer 
arithmetic  cause  inaccuracies.  It  should  be  noted,  however,  that 
the  series  method  is  a  straightforeward  and  quite  practical  method 
for  computing  ^  provided  the  moduli  of  the  numerator  of  the  first 
term  of  the  series  do  not  exceed  25. 
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The  accuracy  of  the  result  may  be  assured  in  either  one  of  two 
ways 4  Each  term  computed  may  be  examined  apd  whenever  the  absolute 
value  of  each  element  of  the  matrix  of  that  term  is  less  than  the 
required  accuracy  computation  may  be  considered  complete*  A  less 
difficult  method  from  the  standpoint  of  programming  is  the  use  of 
a  procedure  set  forth  in  Section  4-19  of  Elementary  Matrices  (Ref  3:145)« 
The  upper  bound  for  the  power  of  a  matrix  is  first  computed  by  this 
procedure*  The  result  is  then  divided  by  ! »  One  thus  computes 
the  upper  bound  of  the  matrix  for  each  term  of  the  series  and 

/tft 


compares  it  with  the  accuracy  required.  When  the  upper  bound  falls 
below  the  required  accuracy  the  number  of  the  term  just  calculated 
is  read  out  and  employed  to  terminate  the  rest  of  the  programs 
involving  series  computations* 

Computation  of  the  fundamental  matrix  for  a  problem  of  a  practical 
nature  is  usually  not  feasible  by  the  series  solution  method  because 
of  the  large  number  of  terms  necessary  for  accuracy.  However,  another 
method  involving  direct  numerical  integration  of  the  system  equations 
can  be  employed* 

Computation  of  ^  .  g-  .  and  tl  bv  the  Method  of  Runoe-Kutta 
The  matrices  p ,  a,  and  can  be  computed  by  direct  integration 
of  the  system  equations*  It  is  the  purpose  of  this  section  to  \  ' 

develop  the  theory  for  this  method*  Consider  the  solution  to  the 
system  equations: 


rcr)  ^  i(T)^  r.  ->■  ^ir)Ur 


(72) 
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Examination  of  equation  (72)  reveals  that  in  the  linear  system  the 
solution  is  constructed  by  superposition  of  the  response  due  to  the 
initial  condition  X(c)‘  C  and  the  response  due  to  the  input 
Taking  first  of  all 

£=^Lt)C  (95) 

one  notices  that  this  could  be  regarded  as  the  solution 

W)  =  ^cr)c  +  P  (125) 

Equation  (125)  is  identical  to  the  solution  for  the  unexcited  system 
with  an  initial  condition 

-  £  (126) 

The  term  "unexcited"  means  that  the  driving  function  is  zero 
in  equation  (72)«  Equation  (125)  is  thus  the  solution  to  the  homo~ 
geneous  system 

=  f  }  y  C^}  ^  (127) 

at  X,~T»  Direct  integration  of  aquation,  (127)  diill^thsn  detensins 
The  solution  for  the  vector  A  will  be  considered  next.  Taking 


the  integral 
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and  carrying  out  the  integration  as  before,  one  has 


--  f fe '-f-Xj  f- 


(128) 


Now  take  the  integral 


(129) 


(13D) 


and  carry  out  the  integration! 

gtf  ^  =  -f-'e-’'-  i 

The  results  of  equations  (128)  and.  Il30)  are  identical.  .The  original 
integrals:  may  thus  be  equated! 

''0 

Referring  back  to  the  solution  to  the  system  equations,  which  is 
equation  (72)  one  sees  that 

g  =  x^)  z  p  +  ■-<'■) I  (i)  cHr  (^32, 


Equation  (132)  says  that  the  vector  ^  may  be  calculated  by  finding 
the  solution  at  time  i,*Z^o  the  system  equations  which  have  been 
excited  with  a  step  input  of  unity  with  an  initial  condition  HCp)  “  0  • 
Thus  the  system 

ri:t)  =  f  ;  xCc>)^0  (133) 
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may  be  integrated  from  X.  -  0  to ''t  to  determine  A  . 

Finally,  one  can  compute  the  matrix  H  by  noting  that 

-  K^f  (134) 

This  can  be  expressed  alternately  as 

%  Kt-jUz)  =  ^  ('*^r)sr  ■+  g  (135) 

'y*i=  (k-i-x)  k-i 

Thus  the  columns  of  the  H  matrix  can  be  calculated  by  computing 
the  state  vector  of  the  unexcited  system  at  the  control  subintervals 
starting  with  an  initial  condition  %(b)  One  then  must  ivaliuat^e 
.the  following  system' at  t  then.', integrate’  from:ii»O  to  C.j.:  t.  , 

from  t  to  2  tT, ,  etc.. : 

^(^6)  ^  f  -  2r  (136) 

Now  that  formulae  ShicK'  can  be  directly  integrated  to  find,  p^.  d.,  ..and  - 
have  been  developed, a  practical  method  of  integration  must  be  considered. 

Again,,  the  only  practical  method  of  integration  is  an  approximate 
numerical  solution.  Of  the  many  methods  available,  the  method  of  Runge- 
Kutta  is  quite  suitable  for  this  type  of  problem  since  it  affords 
accurate  results  and  since  the  computation  procedures  are  recurrent. 

Also,  it  may  be  extended  quite  easily  from  the  scalar  form  to  the 
matrix  form.  Consider  one  equation  of  a  Stati  Space  system  in  sub¬ 
scripted  form: 

^  U37) 
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Replace  by  a  small  interval  and  apply  the  following 
formulae{Ref  10:  123  )* 
if  f  Lxj  ±)  ;  %{0]  =  C. 

than 

cht  - 
Us.  - 

The  formulae  may  in  like  manner  be  applied  to  the  successive  State 
equations  to  compute  In  matrix  form  this  is  more  com¬ 

pactly  expressed  for  the  particular  problem  of  calculating  ^  t 

-XtO)  ~  C  =  £ 

^3  ~  E  ( ^  f 

E  ( 2L  )  -h  ^ 

^2^  .=  i  C^i  ^  ^6^^ ) 
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The  value  of  ^xjust  calculated  is  then  added  to  Xfo)ary6  the  new= 
value  of  X  is  used  as  an  initial  condition  for  another  round  of 
computations.  The  computation  thus  proceeds  in  an  iterative  fashion 
until  the  particular  value  of  ^  is  reached.  The  value  of  2^  at  this 
point  is  then  read  out  as  the  desired  solution. 

It  is  possible,  using  the  method  of  Runge-Kutta  just  described, 
to  compute  P  ,  At  and  H  by  direct  integration  of  the  system 
equations.  The  utility  of  this  method  is  obvious  in  that  it  does  not 
depend  on  a  series  solution  and  hence  the  accuracy  is  not  limited 
by  the  number  of  terms  involved.  Thus,  a  problem  for  which  the  moduli 
of  the  matrix  £  are  large  may  be  handled  as  easily  as  any  other.  Also, 
when  computing  ^  ,  it  is  only  necessary  to  read  out  the  solution 
at  the  end  of  each  successive  control  subinterval,  rather  than  to 
start  from  the  beginning  for  each  column  as  with  the  series  solution. 
One  must  only  be  careful  to  choose  a  proper  value  for  Too  small 

a  ^;6will  result  in  a  roundoff  error  when  using  a  digital  computer 
because  the  ^Xeomouted  will  be  quite  small.  Conversely,  when  ^  ^ 
is  too  large  the  accuracy  t<ill  suffer  because  the  slope  of  the 
function  will  not  remain  constant  over  the  interval  of 

Calculation  of  p,  dt  and  /V  represents  precomputation  which 
is  accomplished  before  the  main  program,  which  is  the  descent,  .-is 
entered.  This  not  only  breaks  up  the  computation  time,  but  also 
allows  one  to  check  for  reasonpble  results  of  the  precomputation 
before  employing  the  descent  program.. 
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The  Descent 

Wort  that  th*  p’ieeomputatlon  has  been  explaihed,:  it  is  only  necessary 
to  list  the  equations  of  the  descent  in  the  order  in  which  they  are 
programmed.  It  should  be  noted  that  an  initial  value  for  the  driving 
-function  40  must  be  guessed  at  and  a  control  interval  T  selected. 

It  goes  without  saying  that  a  good  initial  guess  for  will  greatly 
speed  the  solution.  The  choice  of  the  control  interval  will  of  course 
depend  on  the  linear  transfer  function  itself,  on  the  limits  of  the 
driving  function,  .on  the  initial  conditions,  and  on  the  R  matrix. 

Choice  of  the  R  matrix,  which  is  a  matrix  of  weighting  constants, 
depends  on  the  amount  of  direct  control  it  is  desired  to  exert  on  each 
state  variable*  Further  discussion  of  the  choice  of  the  matrix  jR  will 
be  postponed  until  the  end  of  this  chapter. 

Using  the  precomputed  values  for  p  and  hf  ,  and  a  selected 
the  following  equations  comprise  the  descent  portion  of  the  solution: 


-kCt) 

=  p  +  H  ^ 

(96) 

V 

=  l^Ll)  R  2£(r) 

(61) 

i 

^  H'  R 

(140) 

(141) 

--  -Vt 

(110) 
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The  operator  must  now  choose  between  the  two  descent  schemes.  The 
equations  for  the  Gradient  Projection  Descent  Scheme  are  presented 
first : 

'  ^  <  SA^<[yA- ^ 

[yi- 

A  new  value  of  '7^  must  now  be  computedt 


'>1,  i  .  -  f  ^ 

p  H‘  jk  a  UJ 

The  above  is  then  limited  io  unity: 

^  {yj,) 

The  change  in  the  driving  function  is  calculated  next: 

Finally 

y^^ytuMx-  -  f" 


(142) 


(116) 


(110) 


(123) 


.  If  the  operator  desires,  the  Comer  Aiming, Descent  .Scheme  may  be  used: 

^A-  <f  O  “fek  7 

IF  n>o 


(119) 
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Again  a  new  value  far  '^must  be  calculated: 

"  Sl'^ 

H  '  I<1  ijyr 

and 


From  which 

z  'Vt 


(143) 


(121) 


(122) 


Finally 

(123) 

The  computation  is  repeated  with  the  new  value  of  The 

process  is  iterative  until  the  following  conditions  are  fulfilled 
for  all  control  subintervals,  whereupon  it  is  stopped  and  the  final 
control  is  read  out: 


1  Y/i  -  1  <  S 

FP/t 

<5*,<  -  e 

<  rx 

FPK 

Iftl  <  £ 

(  1  ^  S 

PePK 

>£ 

The  terms  $  and  £’  must  be  used  in  a  digital  computer  solution 
since  the  control  converges  toward  a  boundary,  but  never  quite 
reaches  it.  Also,  the  gradient  never  quite  becomes  zero,  although 
it  approaches  it  so  closely  that  no  noticeable  improvement  can  be  ' 
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made  in  the  terminal  norm.  Thus,  further  computation  would  be  wasted 
effort. 

The  Matrix  S 

A  brief  discussion  of  the  choice  of  the  matrix  R  is  now  in 
order.  The  function 

V=  R  -^cr)  (61) . 

is  termed  a  "quadratic"  if  R  is  symmetrical.  When  this  is  expanded 
one  has 

V  =  %It)  +  'Tgi.  'X^Ct)  f  •  •  •  +  XXlT)  + 

The^'  terms  are  called  mutual  dependance  terms,  while  the  terms 
are  called  self  dependance  terms.  When  formulating  the  criterion  V  , 
it  is  desirable  to  keep  only  the  self  dependance  terms,  since  one 
is  usually  not  interested  in  controlling  the  product  of  two  state 
variables*  Thus,  the //i;^  are  made  zero.  The then  become  weighting 
constants  which  reflect  the  "utility"  of  controlling  directly  each 
state  variable. The  term  "directly^  means  that  an  attempt  is  made 
to  drive  that  particular  state  variable  to  zero  during  the  control 
interval.  This  becomes  more  clear  when  one  notes  that  for  any  /yii 
chosen  to  be  zero  the  corresponding  "JU  will  not  be  driven  to  zero 
but  will  follow  according  to  the  equation 


-  ^<'-i 
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With  these  points  in  mind  the  engineer  may  choose  the  elements  of 

the  matrix  B  in  order  to  solve  the  particular  problem  he  is  confronted 

with* 

As  in  all  minimization  or  maximization  problems  the  choice  of 
the  criterion  used  influences  the  answer  a  great  deal.  Since  the 
matrix  ^  is  a  part  of  the  criterion  used  for  the  particular  method 
discussed  in  this  report  it  will  certainly  influence  the  answer  a 
great  deal,  as  will  be  seen  in  the  next  chapter. 

Determining  the  Optimum  Control  Interval 

After  the  iteration  process  has  been  completed  and  the  optimum 
driving  function  has  been  calculated  for  u  particular  control  interval, 
the  computed  value  for  the  terminal  norm  V  is  read  out.  The  process 
is  then  repeated  for  a  new  control  interval.  After  this  i?  done  several 
times,  the  values  of  V  are  plotted  against  the  respective  control 
intervals  used  for  the  computations.  The  result  is, that  shown  in 
Figure  7  on  the  next  page.  The  end  of  the  solid  part  of  the  curve 
represents  the  last  value  of  V  plotted.  A  tangent  to  the  curve  is  then 
projected  from  this  point  to  the  T  axis  as  shown.  The  value  of  T  at 
the  intercept  of  the  tangent  with  the  T  axis  can  then  be  used  as 
an  estimate  for  a  new  control  interval  for  the  next  computation. 

The  process  may  then  be  repeated  successively.  Note  that  the  dashed 
portion  of  the  V  curve  is  assymtotic  to  the  T  axis.  Since  the  curve 
never  reaches  zero,'  a  small  percentage  of  the  initial  value  for  V/ 
is  picked.  This  value  is  plotted  on  the  V  axis  and  a  horizontal  line 
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is  drawn  as  shown  in  Figurs  8.  When  the  value  of  V  just  falls  below 
this  line  the  corresponding  "T  is  the  optimum  control  interval  and 
the  driving  function  computed  at  this  value  of  T  is  the  optimum 
driving  function. 

This  completes  the  description  of  the  computer  programs.  The 
next  chapter  of  this  report  will  be  devoted  to  the  analysis  of  the 
Minneapolis-Honeywell  relay  servo’ loop. 
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V,  State  Space  Analvaia  of  the  MinneaDolis-HonBvwell  Relay  Servo  Loop 

In  this  part  of  the  report  the  numerical  solutions  for  the  optimum 
driving  function  of  the  linear  controller  portion  of  the  Minneapolis- 
Honeywell  Adaptive  Control  System  will  be  discussed*  In  Chapter. II 
it  was  stated  that  only  the  relay  servo  loop  of  the  Honeywell  system 
would  be  snalyzsd  in  this  report,  since  the  model  could  easily  be 
analyzed  seperately.  A  block  diagram  of  the  relay  servo  loop  is 
presented  on  the  next  page  in  Figure  9* 

,A  meaningful  method  of  analysis  for  any  servo  system  is  to 
find  the  response  to  a  unit  step  input;  this  method  will  be  employed 
in  this  report.  In  Chapter  III  a  method  for  finding  the  optimum  driving 
function  to  drive  an  initial  error  to  zero  in  minimum  time  was 
developed  for  a  linear  controller*  If  the  relay  servo  loop  of  Figure 
9  is  excited  by  a  unit  step  input,  it  will  appear  at  the  output  of 
the  summer  as  a  unit  step  error,  provided  the  system  was  in  equilibrium 
to  start  with*  The  switching  logic  and  the  relay  will  then  generate 
a  driving  function  for  the  linear  controller  to  drive  the  error  back 
to  zero.  Equivalent  to  a  unit  step  input  with  an  initial  output 
error  of  zero  would  be  a  zero  input  with  an  initial  output  error 
of  a  unit  step*  The  output  of  the  summer  would  be  the  same  in  both 
cases  disregarding  the  signs.  Since  the  two  are  equivalent,  the 
problem  that  will  be  considered  in  this  chapter  is  that  of  computing 
the  optimum  driving  functioil  for  the  linear  controller  portion  of  . 
the  Honeywell  relay  servo  loop  when~the  system  has  an  initial  error 
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F igure  9 

Ralay  Servo  Loop  Used  in  Analysis 
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of  a  unit  step* 

The  purpose  of  the  analysis  of  this  chapter  is  twofold:  first, 
to  obtain  an  optimum  driving  function  for  the  linear  portion  of  the 
controller  which  can  ba  compared  to  the  relay  output  of  the  Honeywell 
relay  servo  loop  as;  observed  during  the  analogue  computer  simulation; 

S 

second,  to  demonstrate  the  capabilities  and  limitations  of  the  StSte 
Space  techniques  and  ttka  method  of  Staepest  Descent* 

Solutions  were  attempted  to  five  problems.  Two  of  these  problems 

involved  a  fourth  o'rder  approximation  of  the  sixth  order  linear 

j 

controller.  The  five  solutions  sought  were:  1)  a  solution  for  direct 
control  of  all  six  state  variables  of  the. sixth  order  system,  2) 
a  solution  for  direct  control  of  the  first  four  state  variables 
of  the  sixth  order  system,  3)  a  solution  to  control  dirSctly  all 
four  state  variables  of  the  fourth  order  approximation,  4)  a  solution 
to  control  directly  only  the  first  two  state  variables  of  the  sixth 
order  system,  and  5)  a  solution  to  control  directly  only  the  first 
two  state  variables'  of  the  fourth  order  approximation.  The  fourth 
order  approximation  was  employed  after  numerical  analysis  of  the 
sixth  order  system  indicated  that  the  descent  would  probably  be 
slow  to  converge,  a  fear  that  was  later  borne  out. 

The  Fourth  Order  Approximation 

The  transfer  function  for  the  sixth  order  linear-  controller  is 


(s^f  STAS  4  37')(s‘+Y.rjs  fj.o; 


(147) 
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Examination  of  equation  (147)  indicates  that  a  good  approximation 
to  the  system  could  be  made  by  eliminating  the  high  frequency  poles 
and  zero.  The  result  is  the  transfer  function 

.  -7.1  (S-h.JS)  L^^-S) _  (140) 

Also,  compare  Figure  10  with  Figure  11.  The  former  is  the  Log  Magnitude 
and  Phase  Angle  Diagram  of  the  sixth  order  transfer  function,  while 
the  latter  is  the  Log  Magnitude  and.  Phase  Angle  Diagram  for  the 
fourth  order  transfer  function.  Note  that  the  phase  margin  of  45° 
occurs  at  about  the  same  frequency  on  both  diagrams  as  does  the 
gain  crossover.  The  fourth  order  transfer  function  thus  represents 
a  reasonable  approximation  of  the  Honeywell  linear  controller  at 
the  frequencies  which  will  be  pertinent  to  this  analysis. 

The  Minneapolis-Honevwell  Switching  Logic 

In  this  section  the  Honeywell  switching  logic  will  be  analyzed 
in  order  to  determine  a  set  of  weighting  constants  to  be  used  in  the 
matrix  R  «  The  transfer  function  of  the  switching  logic  is 

e,  _  >^-22S  s 

e  ^  I  +  o/s  5 

If  this  transfer  function  is  excited  by  a  unit  ramp  input  the 
output  is  essentially  proportional  plus  derivative.  For 
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the  output  in  the  time  domain  becomes  approximately 


(151) 

When  the  transfer  function  for  the  linear  controller  portion 
of  the  Honeywell  relay  servo  loop  is  expressed  4-n  State  Space  notation 
where  %i  is  the  system  error.  The  switching  logic  then 
forms  the  sum  .  when  Tt,  alone  is  applied  to  the  switching 

logic.  Note  that  when. driving  the  controller  using  the  Honeywell 
switching  logic  only'  the  first  two  state  variables  are  measured  and 
used  as  a  basis  for  direct,  control.  The  other  four  state  variables 
are  not  directly  controlled,  but  follow  according  to  the  aquation 


^±r.L  4  ^ 


(146) 


Finally,  viewing  the  switching  logic  from  the  standpoint  of  formu¬ 
lating  a  terminal  norm,  one  could  say  that  weights  assigned  to  X.( 
and  T^jfor  direct  control  are  1  and  ,225  respectively. 


Choice  of  the  Matrix  R 

In  Chapter  IV  it  was  stated  that  the  choice  of  the  matrix  J3 
depended  on  the  "utility"  of  controlling  the  various  state  variables 
directly.  Admittedly  the  Honeywell  switching  logic  did  not  employ 
the  quadratic  criterion  function,  but  it  rather  summed  algebraicly 
the  two  weighted  state  variables.  Nevertheless,  to  form  a  basis 
for  comparison,  the  //f,  and  the of  the  matrix  were  chosen 
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to  be  5  and  1.12  respectively.  The  rest  of  the  ^  matrix  was  chosen 
to  be  zero.  Two  problems  were  solved  using  this  matrix  for  R  :  one 
for  the  sixth  order  transfer  function,  and  one  for  the  fourth  order 
approximation. 

In  order  to  demonstrate  the  limitations  and  capabilities  of  the 
methods  used,  three  other  problems  were  formulated.  Two  solutions 
to  the  sixth  order  problem  with  S  diagonal  were  attempted.  In  the 
first  case  the  diagonal  elements  were  6,.  5,  4,  3,  2,  and  1  starting 
with  ;  in  the  second  case  they  were  4,  3,  2,  1,  0,  and  D.  Another 
solution  which  was  attempted  was  one  to  the  fourth  order  approximation 
with  diagonal;  the  diagonal  elements  were  4,  3,  2,  and  1. 

Data  for  the  Solution 

Th^  sixth  and  fourth  order  transfer  functions  were  put  into 
State  Space  form  by  means  of  the  linear  transformations  developed 
in  Chapter  III.  Tables  I  and  II  on  the  following  two  pages  show  the 
data  used  in  the  computer  solutions  for  the  sixth  and  fourth  order 
problems  respectively.  The  number  of  control  subintervals,  and  hence 
the  lengths  of  the  control  intervals  used  were  varied. 

Note  the  magnitude  of  the  last  two  elements  of  the  ^  vector  in 
the  sixth  order  problem.  This  is  due  to  the  nature  of  the  recurrence 
formula 
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Table  II 

Data  for  Fourth  Order  Problem 
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when  applied  to  eliminate  right  hand  side  derivatives.  For  the  sixth 
order  transfer  function  the  coefficients  of  the  denominator  were 
quite  large  which  led  to  a  rapid  buildup  in  magnitude  of  the  elements 
of  ^  ,  Note  that  the  ^  vector  for  the  fourth  order  problem  does  not 
exhibit  this  characteristic.  Examination  of  the  matrix  for  the 
sixth  order  problem  (Table  Oil,  Appendix  D)  shows  that  the  large 
values  of  the  last  two  elements  of  ^  lead  to  large  values  for  part 
of  the  last  two  columns  of  Experience  with  computation  shows 
that  these  large  variations  in  the  elements  of  ^  cause  a  rather 
slow  descent.  The  reason  for  this  can  be  seen  from  the  formula 

'V)  r  ,  (141) 

fU'  S 

The  large  values  of  H  relative  to  ^  cause  'Iff  to  become  quite  small 
and  thus  the  change  in  control  becomes  very  small.  This  last 
fact  also  holds  for  other  calculations  of  '»7» 

Results  of  the  Computation 

In  this  section  the  results  of  the  computation  will  be  discussed 
All  of  the  figures  referred  to  in  this  section  are  located  at  the 
end  of  this  chapter#' 

Figure  12  is  a  plot  of  the  optimum  driving  function  for  the 
fourth  order  approximation  with  all  four  state  variables  directly 
controlled.  Figure  13  presents  a  plot  of  the  state  variable  trajec¬ 
tories  which  result  when  the  control  of  Figure  12  is  applied.  A 
plot  of  the  terminal  norm  V  for  the  different  control  intervals  used 
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is  presented  in  Figure  14.  The  optimum  control  interval  for  a 
terminal  norm  V  of  .082  is  1,85  seconds*  Three  results  should  be 
pointed  out.  First,  the  control  is  essentially  "bang-bang”.  Secondly, 
the  state  variables  are  not  all  derivatives  of  one  another.  This  can 
be  immediately  explained  by  noting  that 


Hi  ^ 

but  that 


(152) 


(146) 


Thirdly,  note  that  the  negative  excursion  of  is  slightly  greater 
in  magnitude  than  the  original  error, 

A  solution  for  the  sixth  order  system  with  weighting  constants 
of  6,  5,  4,  3,  2,  and  1  was  not  completed.  The  descent  was  attempted 
but  after  30,000  iterations  data  showed  that  convergence  was  going 
to  be  quite  slow.  There  is  some  question  as  to  whether  control  is 
possible  at  all  for  .this  case.  The  reasons  for  this  will  become 
more  clear  after  the  data  from  two  of  the  other  problems  is  reviewed. 
The  data  from  the  solution  for  direct  control  of  the  first  two 
state  variables  for  both  the  fourth  order  and  the  sixth  order  systems 
is  to  be  considered  next.  The  various  plots  of  the  optimum  control, 
of  the  state  variables,  and  of  the  terminal  norm  are  shown  in  Figures 
15  through  24.  Note  the  differences  in  scale  which  were  necessary 
for  the  plots  of  the  sixth  order  problem.  The  control  interval 
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for  the  V  of  the  fourth  order  approximation  to  reach  a  value  of  .16 
is  .65  seconds.  In  .675  seconds  the  terminal  norm  of  the  sixth  order 
system  reached  a  value  of  .13.  The  relay  switching  occurs  for  both 
systems  at  about  the  same  time.  Note  also  the  close  resemblance  of 
the  trajectories  of  and  Tt^for  both  systems.  These  resemblances 
suggest  that  the  fourth  order  system  was  a  good  approximation  of 
the  sixth  order  system.  Of  special  interest,  however,  are  the  tra- 
Jectoxes  of  the  higher  order  state  variables.  Remember  that  direct 
control  was  not  attempted  for  the  state  variables  through  TCft 
The  plots  of  the  state  variable  trajectories  for  the  sixth  order 
system  illustrate  quite  well  what  has  happened.  The  State  Space 
equations  for  the  sixth  order  system  show  that 


but  more  especially  that 

IT' 


- 

-  — - 
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(154) 


(155) 
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The  trajectories  of  the  state  variables  show  these  relationships 
graphically.  Note  sspecislly  how  the  trajectories  of  7^, 

“Xu  governed  quite  apparent  why  direct  control 

of  all  the  state  variables  is  so  difficult  to  obtain  in  this  problem* 
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The  control  is  fighting  itself  in  this  case.  It  should  be  pointed 
out  again  that  the  presence  and  size  of  the  s  is  a  result  of 
the  numerator  dynamics  of  the  entire  system  and  the  high  natural 
frequency  of  the  servo  and  actuator. 

The  trajectories  of  and  'J^for  the  fourth  order  example 
cannot  be  compared  with  those  of  Ti^and  Xf.for  the  sixth  order  system, 
again  because  of  the  different  way  in  which  the  numerator  dynamics 
influence  the  solution. 

A  solution  for  direct  control  of  the  first  four  state  variables 
of  the  sixth  order  system  was  also  attempted.  Convergence  of  the 
solution  was  quite  slow  and  therefore  a  complete  solution  was  not 
obtained.  However,  it  appears  that  convergence  would  be  possible 
if  enough  iterations  were  tried.’ 

A  term  which  seems  applicable  to  a  problem  of  this  nature  is 
"controllability".  A  definition  of  "controllability"  is  not  easily 
stated,  but  it  is  possible  to  think  of  it  in  terms  of  the  number 
of  iterations  in  the  descent  necessary  to  compute  the  optimum  control. 
The  more  iterations  necessary  the  less  controllable  is  the  system. 

Of  course  it  is  immediately  evident  that  the  controllability  depends 
on  which  of  the  state  variables  one  wishes  to  control  directly  and 
to  what  degree.  For  example,  controlling  directly  all  six  state 
variables  of  the  sixth  order’ system  to  the  same  degree,  would  be 
almost  impossible.  Conversely,  the  descent  for  direct  control  of  the 
first  two  state  variables  of  the  sixth  order  system  was  completed 


GE/EE/62-8 


quite  rapidly.  See  Tables  D-III  through  &-V  in.  Appendix  0  for  a 
comparison  of  the  number  of  iterations  necessary  for  the  different 
problems. 

The  results  of  this  section  have  demonstrated  the  use  of  the 
State  Space  techniques  and  the  Steepest  Descent  method  in  solving 
a  practical  problem  of  computing  an  optimum  driving  function.  An 
attempt  was  made  to  demonstrate  the  limitations  of  the  methods  ae 
well  as  the  capabilities.  Data  which  can  be  correlated  with  the 
analogue  computer  simulation  of  the  problem  was  presented.  Attempts 
at  five  solutions  to  the  problem  were  made;  three  of  these  were 
successful. 

The. next  part  of  this  report  will  be  concerned  with  the  analogue 
computer  simulation  of  the  Minneapolis-'Honeywell  relay  servo  loop. 
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VI.  Analogue  Computer  Simulation  of  the  Minneapolis-Honeywell 

Relay  Servo  Loop 

In  order  to  form  a  basis  for  comparison  of  the  results  of 
Chapter  V  an  analogue  computer  simulation  of  the  Minneapolis-Honey— 
well  relay  servo  loop  was  performed.  Figure  9,  which  is  a  block 
diagram  of  the  relay  servo  loop  is  presented  at  the  end  of  this 
chapter  for  the  reader's  convenience. 

In  keeping  with  the  analysis  of  Chapter  V.  the  model  was  eliminated 
from  the  original  Honeywell  Adaptive  Control  System,  as  were  the 
gain  changer  and  limiter  and  the  AC  dither.  The  response  of  the 
relay  servo  loop  with  a  step  input  applied  is  that  which  is  desired. 

The  relay  output  and  the  system  output  can  then  be  compared  with 
optimum  driving  function  and  the  state  variable  trajectories  obtained 
by  the  calculations  of  Chapter  V« 

Description  of  the  Simulation 

The  simulation  of  the  relay  servo  loop  is  quite  straightforeward 
arid  will  not  be  discussed  in  detail. The  procedures  outlined  in 
Chapter  19  of  Aotomatic  Control  iSvs'^feros  Analysis  and  Synthesis  - 
were  followed. (Ref  2:453-475). 

The  following  points  concerning  the  simulation  should  be  noted 
howeveri  l)  due  to  the  high  natural  frequencies  of  the  servo  and 
acutator  and  the  rate  gyro  it  was  necessary  to  time  scalattha  problem 
by  slowing  the  solution  by  a  factor  of  ten|  2) the  electronic  relay 
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used  in  the  Honeywell  loop  was  replaced  by  a  high  gain  amplifier 
in  cascade  with  a  differential  relay.  Figure  26  at  the  end  of  this 
chapter  is  a  diagram  of  the  simulation, 

A  describing  function  analysis  assuming  an  ideal  relay  with 
an  output  of  1,5  volts  was  used  to  determine  stability  and  the 
approximate  frequency  of  the  limit  cycle  to  be  expected,  A  Log 
Magnitude  Phase  Angle  Diagram,  for;  the  describing  function  analysis 
is  shown  in  Figure  27  at  the  end  of  the  chapter. 

Results  of  the  Simulation 

Step  inputs  of  one  volt  were  applied  to  the  system  after  a 
steady  state  limit  cycle  had  been  established.  The  one  volt  corresponds 
to  the  one  unit  of  initial  error  used  in  the  computation  of  the 
optimum  driving  function  in  Chapter  V.  Figure  .25  at  the  end  of  this 
chapter  is  a  typical  analogue  computer  recording  for  the  response 
to  the  unit  step  input.  The  traces  are,  beginning  at  the  top,  error 
signal,  input  to  the  high  gain  amplifier,  relay  output,  and  system 
output,  in  that  order.  Upon  application  of  the  step  input,  initial 
switching  of  the  relay  took  place,  if  the  polarity  wasn’t  already 
such  as  to  drive  the  error  to  zero.  As  the  relay  input  crossed  the 
aero  axis  the  relay  switched  polarity.  Finally,  when  the  error 
approached  zero  limit  cycle  operation  was  resumed. 

Results  of  high  speed  recorder  runs  show  that  the  error  is 
reduced  to  zero  in  about  ,75  seconds.  The  limit  cycle  is  resumed 
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after  about  .65  seconds.  The  average  time  from  the  detection  of  the 
initial  error  until  relay  pelarity  reversal  is  bout  .4  seconds. 

Exact  times  are  difficult  to  determine  because  they  depend  upon 
the  state  of  the  system  in  the  limit  cycle  to  some  extent. 

It  should  again  be  pointed  out  that  the  input  to  the  high  gain 
amplifier  is  the  error  signal  plus  .225  of  its  derivative.  The  trace 
of  the  error  signal  shows  that  for  a  negative  error  the  slope  is 
positive;  thus,  when  the  error  decreases  below  .225  of  its  derivative 
at  that  point,  the  relay  reverses  polarity.  This  reversal  prevents 
overshoot.  Finally,  note  that  the  shape  of  the  trace  of, the  system 
error  resembles  the  shape  of  the  trajectory  of  "X.*  in  Figure  16  which 
is  located  in  the  preceding  chapter. 

The  data  just  presented  leaves  little  more  to  be  said  concerning 
the  analogue  computer  simulation.  The  reader  should,  however,  keep 
the  results  in  mind  since  they  will  be  compared  in  the  next  chapter 
of  this  report  with  the  results  of  the  computations  for  optimum 
control. 
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VII.  Conclusions  and  Recommendations  ■ 

In  this  last  chapter  of  this  report  the  data  compiled  in  the 
previous  sections  will  be  summarized  and  conclusions  will  be  drawn 
from  it.  Recommendations  for  further  study  of  certain  phases  of 
this  report  will  then  be  made. 

It  is  not  the  purpose  of  this  chapter,  or  of  this  report,  to 
attempt  to  determine  the  quality  of  the  Minneapolis-Honeywell  Adaptive 
Control  System  just  discussed.  Rather  it  is  to  demonstrate  the  ap¬ 
plication  of  State  Space  and  optimization  techniques  by  analyzing 
a  practical  problem. 

Summary  of  Data 

In  this  section  data  from  the  previous  chapters  will  be  summarized. 
Also,  from  time  to  time,  the  figures  of  Chapters  V  and  VI  will  be 
referred  to  without  stating  their  location. 

It  is  possible  to  compare  only  part  of  the  data  from  the  computations 
for  the  optimum  driving  functions  to  that  from  the  analogue ' computer 
simulation.  The  data  that  may  be  compared  are  the  results  of  the 
optimum  solution  to  the  fourth  and  sixth  order  problems  for  direct 
control  of  the  first  two  state  variables.  A  glance  at  Figures  15 
and  22  shows  that  the  driving  function  calculated  compares  very  well 
with  the  output  of  the  relay  as  recorded  during  analogue  simulation. 

The  first  relay  polarity  reversal  for  both  computational  solutions 
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occured  at  about  .4  seconds,  as  did  the  first  polarity  reversal 
during  analogue  simulation.  The  simulation  also  showed  the  error 
driven  to  essentially  zero  in  about  .7  seconds.  This  corresponds- 
quite  well  to  a  terminal  norm  of  ,13  with  a  control  interval  of 
.675  seconds.  Of  course  it  is  impossible  to  state  that  a  direct 
correlation  exists  because  the  criterion  formed  by  the  switching 
logic  of  the  Honeywell  relay  servo  loop  and  that  formed  by  the 
terminal  norm  function  are  different.  However,  the  simularity  of 
the  two  results  does  indicate  that  for  this  particular  problem  the 
criterion  function  used  for  the  computation  of  the  optimum  driving 
function  is  meaningful. 

Note  the  shape  of  the  trajectory  of  'JC/  which  is  plotted  in 
Figure  16.  It  is  about  the  same  as  that  of  the  trace  of  the  error 
signal  which  was  recorded  during  the  analogue  simulation. 

Another  point  of  interest  in  the  results  of  the  digital  com¬ 
putation  is  the  shapes  of  the  trajectories  of  the  various  state 
variables.  Comparison  of  Figure  13  with  Figure  23  shows  that  when 
an  attempt  was  made  to  directly  control  and  ^in  the  fourth 
order  approximation  their  'excursions  were  greater  during  the  control 
interval  than  when  they  were  allowed  to  follow  without  direct  control 
Of  course  X^andX»were  no.f  driven  to.zefcp  in.  the  .latter,  case. 

Another  point  of  interest  is  that  a  longer  control  interval 
was  necessary  to  control  directly  four  state  variables  than  to 
control  two;  thus,  the  expenditure  of  energy  by  the  system  had  to 
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be  greater  for  the  longer  control  interval.  Also,  one  sees  from 
Figure  13  that  the  maximum  negative  excursion  of  was  greater 
than  the  initial  error  in  the  system. 

A  brief  glance  at  the  tabular  data  in  Appendix  D  shows  the 
number  of  iterations  necessary  to  attain  convergence  varies  over 
a  wide  range  for  the  different  problems  considered.  Generally,  but 
not  in  all  cases,  the  number  of  iterations  necessary  increases  as 
the  optimum  control  interval  is  approached.  Also,  the  more  difficult 
it  is. to  control  certain  state  variables,  the  greater  is  the  number 
of  iterations  necessary  for  convergence.  In  general,  then,  the  rate 
of  convergence  is  a  measure  of  the  "controllability”  of  the  system'. 
The  rate  of  convergence  may  be  measured  by  the  number  of  iterations 
necessary  to  attain  convergence. 

Conclusions 

On  the  basis* of  the  data  presented  in  the  preceding  section 
and  that  presented  in  discussions  elsewhere  in  this  report  several 
interesting  conclusions  may  be  drawn.  Probably  the  most  important 
and  far  reaching  of  these  is  that  each  solution  for  the  optimum 
driving  function  is  in  reality  a  sub-optimization.  The  results  of 
Chapter  V  showed  how  much  the  optimum  control  interval  and  the 
trajectories  of  the  state  variables  depended  on  the  matrix  R  .  It 
goes  without  saying  that  a  different  function  for  the  'terminal  norm 
would  also  have  influenced  the  results  a  great  deal.  It  should  be 
understood,  however,  that  the  computed  driving  function  is  optimum 
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for  the  terminal  norm  chosen.  The  question  remains  of  whether  the 
terminal  norm  used  is  the  beat  criterion. .Thus  the  problem  of 
picking  a  criterion  to  minimize  still  remains  a  matter  of  judgement. 
For  example,  consider  the  particular  system  analyzed  in  this  paper. 
Compare  the  solution  for  the  fourth  order  system  with  the  diagonal 
of  R  equal  5,  1.12^  0,  □,  and  that  for  the  same  system  with  the 
diagonal  of  R  equal  4,  3,  2,  1.  If  one  is  interested  in  reducing 
the  error,  which  in  this  case  is  the  first  state  variable,  .to  zero, 
it  appears  better  to  attempt  to  control  directly  only  the  error 
and  its  derivative  than  to  attempt  to  control  directly  all  the 
state  variables.  Certainly  the  optimum  control  interval  is  less  for 
the  former. 

It  should  again  be  emphasized  that  the  criterion  function  chosen 
is  reasonable,  although  some  of  the  values  used  for  the  ^  matrix 
might  have  been  a  bit  unrealistic.  A  further  point  along  these 
lines  is  that  attempts  to  control  state  variables  which  are  dominated 
by  i^Awill  probably  meet  with  very  slow  convergence  or  even  lack  of 
convergence. 

The  question  of  "controllability”  and  the  related  computation  tme 
is,  for  this  report,  somewhat  academic^  except  in  cases  where  it  ia 
impossible  to  obtain  convergence  without  many  thousands  of  iterations. 
However,  one  problem'  for  which  computation  time  would  be  important 
is  that  of  an  adaptive  state  vector  control  system.  Basically,  a 
system  of  this  type  would  operate  in  the  following  way.  The  state 
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variables  wauld  continuously  be  measured.  Consider  that  they  represent 
the  error  of  the  system.  The  necessary  input  to  a  linear  controller 
to  drive  these  variables  to  zero  would  then  be  continuously  computed 
in  real  time  and  applied  to  the  controller.  For  a  linbar  controller 
with  good  controllability  such  computation  in  real  time  seems  feasible 
by  the  Steepest  Descent  method  provided  a  very  high  speed  computer 
is  employed. 

Finally,  the  fact  that  all  the  solutions  for  the  optimum  driving 
function  for  the  examples  investigated  in  this  report  were  "bang- 
bang"  is  of  interest.  The  same  results  were  obtained  whether  the 
Gradient  Projection  Descent  Scheme  or  the  Corner  Aiming  Descent 
Scheme  was  employed.  This  result  is  in  direct  support  of  the  maximum 
principal  of  Pontriaginf  in  addition,  it  should  be  pointed  out  that 
the  "bang-bang"  solutions  were  arrived  at  by  an  entirely  different 
mathematical  approach  than  was  the  proof  of  the  maximum  principal. 

Recommendations 

Unfortunately,  there  were  many  interesting  phases  of  the  problem 
left  untouched  in  this  investigation  due  to  the  element  of  time. 

Also  it  should  be  pointed  out  at  this  time  that  the  computer 
programs  used  and  a  brief  explanation  of  each  are  presented  in 
Appendix  A.  Hopefully  enough  information  has  been  provided  for  the 
interested  reader  to  utilize  the  programs  with  little  added  study. 

•  One  of  the  phases  of  the  analysis  which  remains  unexplored  is 
further  experimentation  with  the  matrix  R  .  Along  this  same  line 
is  the  testing  of  different  criterion  functions  used  in  the  descent. 
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One  criterion  function  of  particular  interest’ would  be 

r  ^yc:,jc^cr)  (157) 

57/ 

The  computation  of  the  optimum  driving  function  to  be  applied 
to  a  non-linear  controller  is  also  possible  by  the  method  of  Steepest 
Oesceht*  The  interested  reader  is  referred  to  Part  V  of  Reference  5« 

In  the  field  of  digital  computer  programming  it  would  be  quite 
useful  to  modify  the  existing  descent  program  to  seek  out  the  optimum 
control  interval  automatically.  Of  course  due;.to  the  computation  times 
involved  the  use  of. such  a  program  would  have  to  be  confined  to  systems 
with  a  reasonable  amount  of  controllability. 

Finally,  a  study  to  determine  the  feasibility  of  using  the 
Steepest  Descent  method  in  an  adaptive  state  vector  control  system 
could'  be  investigated. 

The  above  list  of  possible  topics  in  connection  with  this  report 
is  by  no  means  complete.  However,  it  is  hoped  that  some  stimulating 
ideas  were  presented. 

In  this  report  State  Space  techniques  and  the  Steepest  Descent 
method  were  employed  to  obtain  the  optimum  driving  function  for  a 
practical  linear  controller  in  a  relay  servo  loop.  It  is  felt  that 
the  results  obtained  demonstrate  the  limitations  and  capabilities 
of  this  type  of  analysis  quite  clearly. 
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Appendix  A 

Lietlno  and  Explenatlon  of_the  Dioitel  Cgwputer  Proorewe 

The  puxpoee  of  this  ^pendlx  la  to  preaent  lietlnge  of  the  digital 
computer  pxograma  uaad  for  the  computation  of  the  optimum  control  to> 
gather  with  a  brief  explanation  of  each  program*  It  ia  aaeumed  that 
the  reader  has  a  knowledge  of  the  FORTRAN  progrmeming  ayatem.  In  writing 
the  programs  several  changes  in  notation  wars  nacssaary.  Table  A-I  on 
the  next  page  ia  a  list  of  each  important  symbol  used  in  the  programming 
together  with  the  correaponding  syaibol  used  in  Part  IV.  The  reader  may 
wish  io  refer  to  this  table  from  time  to  time  in  ease  of  confusion* 

Calculation  of  ALIT  bv  Runqa-Kutta 

In  the  first  program  which  ia  listed  on  pagas  164  and  103  ^  is 
calculated  by  the  method  of  Runge»Kutta*  This  method  is  developed  in 
Part  IV  of  this  report*  The  program  is  written  in  1620  AFIT  FORTRAN. 

The  following  explanation  refers  to  the  Hating  of  the  program* 

The  values  for  Att  for  the  order  of  the  system,  and  for  the  control 
subinterval  are  first  read  in  as  DELT,  N,  and  TAU  rsspsctively*  The  p 
and  ^  matrices  are  then  read  in  by  columns*  Statements  5  through  14 
perform  the  Runge-Kutta  calculation  for  A  teat  is  then  mads  to 
aee  if  the  control  aubintarval  TAU  has  been  reached*  If  it  has  not, 
control  ia  ratumed  to  statement  5  where  the  calculation-  is  performed 
again  with  the  new  value  of  TU  Whan  TAU  is  finally  rsachad  g.  is  read 
out  and  punchad  on  carda* 
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Table  A- I 

Meanings  for  Computer  Symbols. 


Computer  Symbol 

Symbol  in  Text 

N 

BELT 

A  yt 

TAU 

Z 

F 

F 

G 

t 

DX 

A  (Descent  Program 

J?  fr  FCA-c) 

A  (Program  for  H ) 

Ou  /»  TClo) 

J£ 

KA 

K-l 

H 

a 

U 

UHI 

ULO 

Yjf 

EPSIL 

a 

DELTA 

€ 

C 

4 

P 

R 

XT 

1L(T) 

W 

ETA  or  ETA2 

'V,  'fl,  gf' 

V 

V 
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CALCULATION  OF  'ALT T . BY"  RONXsr'  ICOTTA' - 

READ»DELT»N*TAU 

DIMENSION  F(6»6) ♦G(6) ♦XNO(6) tDXIS.A) 

DO  1  K»1»N  . . .  . . 

DO  1' I«1»N 

1  READ»F(I»IC) 

DO  2  K»1.N  .  ■■■ 

2  XNO(K)«0. 

DO  3  K»1»N 

3  REAO.G(X)  .  .  . . 

P  =  0. 

5  DO  7  I«1»N 

z*o«  .  '  '  ‘  . . 

DO  4  K«1»N 

4  Z=Z+F(I»K)*XNO(K)*DELT 

K*I  . . .  . .  . . 

7  OX(K»l  )«Z+G(IC)*DELT 

DO  9  I=«1»N  _ 

DO  8  K»1*N 

8  Z»Z+F(  I  .lO^IXNOIXJ+OXlKtD/Z.I^DELT 

K-I  .  .  . . .  ■  . . .  . . . 

9  DX(IC*2)*Z+6(K)*DELT 
DO  11  I«1»N 

Z-0,  - - - 

DO  10  K»1»N 

10  Z«Z+F(  I  ♦K)*{XNO(IC)+OX(K.2)/2.)*DELT 

K«I  ■  . .  . . . 

11  DX(K,3)-Z+G(X)*OELT 
DO  13  I-l.N 

■Z*0.  . — - - - - - 

DO  12  K-l.N 

12  Z«Z+F(  I  .»O*(XN0(X)+DX(Kf3)  )*DELT 

K=I  . .  ■  . . .  ■■■ 

13  DX(K,4)»Z+6(K)*DELT 
DO  14  K«1»N 

14  xHo  (X )  «xNO{ x)+ riDx  (x:.Tr+7;wx-rr»Ti +T;"xoxTTrr3iTi^^ 

p«i.+p 

IF(TAU“P*DELT)  15,15.5 

15- TYPE  16  •  - '  . . 

16  FORMAT  I/18H  FOLLOWING  IS  ALIT) 

DO  17  K«1»N 

. PUNCH  45»XNO('XJ -  "  ~ 

17  TYPE  44.*XNO(IO 
44  FORMAT  (/E14.8) 
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45  FORMAT  (2X,E14.8) 
STOP 
END 
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Caleulatlon  of  Large  A  and  H  bv  Runae-Kutta 

The  program  for  the  calculation  of  £  and  ^  by  the  Runge-Kutta 
method  is  listed  beginning  on  the  next  page.  There  are  two  forma  of 
the  program  included •  one  written  in  1620  AT IT  FORTRAN  and  the  other 
written  in  7090  FORTRAN.  The  only  differences  in  the  two  are  the  input 
and  output  statements  and  the  method  of  controlling  the  use  of  the 
options.  VAien  it  is  desired  to  use  a  large  number  of  control  subintervals 
and  the  order  of  the  problem  is  larga  it  is  advantagsous  to  use  the 
7090  version  aincs  computation  time  is  substantially  less.  Ths  7090 
varaion  will  be  discussed  in  this  appendix. 

The  program  contains  two  options  which  axe  controlled  by  entering 
a  value  of  1  or  0  for  KSENS.  These  options  axe  the  computation  of  H 
or  the  computation  of  jp  (4t)«  It  is  necessary  to  calculate  the  latter 
for  each  control  subintervid  fox  use  in  a  later  progran  which  ealculatea 
the  trajectories  of  the  state  variables. In  case  the  first  column  of  F 
is  zero,  which  corresponds  to  an  integration,  in  tha  linear  controller,  p 
will  remain  constant  throughout  the  control  interval;  thus  its  calculation 
is  unnecessary*  For  this  xaason  the  calculation  of  ^  was  included  in 
this  program  only  as  an  option. 

The  problem,  but  not  the  program,  is  also  restartable.  It  is  only 
necessary  to  enter  tha  last  computed  column  of  ^  as  A(K)  and  start 
ths  program.  As  many  more  columns  of  ^  as  are  desired  will  then  be 
computed.  Finally,  whan  using  the  Runga-Kutta  method  one  usually  starts 
at  -t-  O ;  thus  the  computation  of  the  columns  of  H  starts  at  the  right 
and  proceeds  to  the  left. 
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C  CALCWtATt'OK  -Of  UARtSt  A  ANTT  H'  BY  ROIfOf  TUTTA 

READ  INPUT  TAPE  2 ♦ 99 ♦ N » JE ♦ KA .TAU tDELT »T 
READ  INPUT  TAPE  2.122.K.SENS 
122  FORMAT! 18) 

99  FORMAT(3I4.3F8.0) 

DIMENSION  H(6»99) »F(6.6) fAIS) .DX(6*4) 
TAUl-TAU 
DO  1  K*1»N 
DO  1  I»1*N 

1  READ  INPUT  TAPE  2»112»F(I»K) 

111  FORMAT{E16.8) 

112  F0RMAT(F11.0» 

DO  2  K»1*N 

2  READ  INPUT  TAPE  2»111.A(K) 

P  =  0. 

IF(<SENS)61.62 

61  PEaO. 

63  FORMAT(21H  FOLLOWING  IS  LARGE  A) 

WRITE  OUTPUT  TAPE  3»63 

GO  TO  31 

62  PE»1. 

WRITE  OUTPUT  TAPE  3 *64 

64  FORMAT {25H  FOLLOWING  IS  H  TRANSPOSE) 

31  DO  3  K*1 .N 

I*K 

'  3  H(  I*JE)=A(IC) 

4  DO  6  I«1*N 
Z  =  0. 

DO  5  K»1*N 

5  Z=Z+F( I )*A(K)*DELT 

K=I 

6  DX(l(fl)«Z 
DO  8  I«1*N 

Z  =  0.  . . 

DO  7  k:»i*n 

7  Z»Z+F(  I  .K)*(A()C)+DX(K»1  )/2.  )»DELT 

K»I 

8  DX{K»2)»Z 
DO  10  I«1*N 
Z  =  0. 

DO  9  IC»1*N 

9  Z  =  Z+F{  I  *K  )»(A()C)+DX(IC»2)  /2.  )*DELT, 

)C=I  "  - 

10  0X(K»3)»Z 
DO  12  I»1*N 
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Z-O. 

DO  11  K-l.N 

11  Z«Z+F(  !»<)*(  A  {K»+OX(K*3n*DELT 
'  K«I 

12  0X(<*4)«Z 
DO  13  K-l.N 

13  A  ( K )  -A  ( K )  •*•  ( ( DX  ( K*  1 )  +2;#DXTK  . ZT+Z.^DX-fC. i  f^GXXXt 4T1 76  .  ) 
P«U+P 

131  IF(TAU-P*DELT)  14,l4»4 

14  JE«JE-1  . . . 

PE«1*+PE 

TAU*PE*TAU1 

DO  15  K-1»N  . . . . . . .  . . 

I  »K  - 

15  H( I tJEJ-AIK) 

IF(T-TAU)  16*16»4- . 

16  DO  18  JE-l.KA  • 

18  WRITE  OUTPUT  TAPE  3  ♦  19 »  < HU  ♦  JE )  ♦  I «1  »N ) 

DO  121'  JE«i,KA 

DO  121  I»1»N 

121  WRITE  OUTPUT  TAPE  14 *1 1 1 *8 ( I » JE ) 

19  FORMAT (6E16. 5) 

CALL  EXIT 

END 
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The  order  of  the  syetem,  the  control  eubinterval,  Alt*  the  number 
of  control  eubintervals,  the  control  interval*  and  again  the  number 
of  control  subintervale  are  read  in  as  N*  TAU*  SELT,  JE*  T*  and  KA 
respectively*  F  is  then  read  in  by  statemant  1*  If  H  is  being  calculated 
a  is  read  in  by  atatement  2;  however*  if  p  is  being  calculated*  "Xio) 
is  read  in  by  statement  2*  Statements  £l  and  62  select  the  index 
to  be  used  for  the  columns  of  H  or  p  •  Statament  3  then  begins  the 
Runge-Kutta  calculation  of  AX*  !■  added  to  xfo)  statement  13* 
Statement  131  then  testa  to  determine  whether  the  value  for  TAU  has 
been  reached.  If  it  has  not*  control  is  returned  to  statement  4  and 
the  calculation  of  a  new  value  for^  commences;  if  TAU  has  been  reached* 
the  2^  just  calculated  is  stored  as  a  column  of  y  after  which  the  value 
of  TAU  is  increased*  A  check  is  then  made  to  determine  if  the  new  value 
of  TAU  is  the  same  as  T*. If  it  is  not  control  is  returned  to  statement 
4;  however*  if  TAU  equals  T*  the  matrix  fH  is  punched  on  cards  by  columns 
one  element  to  each  card*  H  is  also  typed  out  as  H^*  Exactly  the  sama 
computation  is  performed  for  p  except  that  different  indexing  and  dif¬ 
ferent  initial  conditions  axe  used* 

Calculation  of  U 

The  coaiputatlon  of  the  optimum  control  is  accomplished  by  the 
program  which  is  listed  beginning  on  page  110. The  program  which  is 

written  in  7090  FORTRAN  contains  two  options.  These  two  are  the  use 
of  the  Comer  Aiming  Descent  Scheme  or  the  Gradient  Projection  Descent 

Scheme.  Refer  now  to  the  program  listing* 
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C  CALCULATION  OF  U 

READ  INPUT  TAPE  2 »3 .N.KA tMAX tUHI »ULO»EPS I L»OELTA 
READ  INPUT  TAP_E  2»5,KSE^NS__ 

DIMENSION  Xt(6)«P<6»6)»Y(6)»bM(6)»A(6) 

DIMENSION  U(99) ♦H(6*99)tC(99) ♦DELU(99) »W<99) .ABSC(99) 
_ DO_2  I»l.ll _ 

2  READ  INPUT  TAPE  2f7,A(I) 

DO  13  JE»1»KA 

13  READ  INPUT  TAPE  2*7,U(JE) 

. 7' FbRMAt(E16.8)  . 

3  FORMAT  ( I4»I4»I7»F8.1.F8.1*F10.7»F7.2) 

_ DO  16  JE»1»KA _ 

DO'  16  I»i»N 

16  READ  INPUT  TAPE  2*7*H(I.JE) 

DO  18  I«1*N _ _ _ 

bo  18  I2»i'»N 

18  READ  INPUT  TAPE  2»4.P(I2»I) 

4  format (F8»4) _ 

ITER»0 

20  DO  24  I«1.N 

_ 

DO  23  JE»i»KA 

23  Z»Z+H(I»JE)#U(JE) 

24  XT(  n»Z4A(I) _ 

V»0* 

'  DO  35  I2«1*N 

Z«0*  _ 

DO  32  i«i»N 
32  Z»Z+P(  I2f  n*XT(  I  ) 

_ UlU _ 

Y(n«z 

35  V»V+XT( I)*z 

_  _ 

DO  41  JE«i»KA 
Z«0, 

_ DO  40  I«1*N _ 

40  Z«Z+H(I.JE)*Y(n 

41  C(JE)»Z 

_ IF(KSENS)^9.199»42 

42  DO  46  JE«1»'KA . 

IF(C(JE))  43.44.45 

43  W( JE)»UHI-U( J EJ _ 

GO'  TO  46 

44  W(JE)»0. 

GO  TO  46 
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45  W(JE)=ULO-U( JE) 

46  D  =  D-C(  JE)*W(  JE). 

47  DO  51  I»1»N 
Z*0. 

DO  50  JE»1.KA 
50  Z»Z+H( I ♦JE)*W( JE) 

5  r”DM  (  I  )  =Z . .  ■ 

Z3«0. 

DO  58  J2=*1*N 
Z2a0. 


DO  55  I«1*N 
55  Z2»22-*-P(  I2»I  )»DM(  1  ) 


58 

I»I2 

Z3»Z3+Z2*DM( I ) 

ETA=D/Z3 

63 

64 

IF(l.-ETA)  66.64*64 

ETA2*ETA 

GO  TO  69 

66 

ETA2«1. 

69 

DO  72  JE«1»KA 

72 

U  ( JE  )  »U  ( JE  )  ■♦•ETA2*W  (  JE  ) 

GO  TO  250 

•  - 

199 

Z=0. 

200 

DO  201  JE-l.KA 

201 

Z»Z-C( JE)*C( JE) 

ETAN«Z 

DO  203  1*1. N 

Z*0. 

DO  202  JE*l.iCA 

202 

Z=Z+H( I .JE)*C(JE) 

203  DM( I )-Z 
ETAD*0. 

DO  205  I2=1»N 


Z*0. 

DO  204  I»1*N 

204  Z«Z-»-P(I2 1 1  )»DM(  I  ) _ _ 

I»i2  “ 

205  ETAO«ETAD+Z*OM( I ) 

ETA*ETAN/ETAD  _ 

53  DO  206  JEViVka  ' 

OELU( JE)»ETA*C( JE) 

W(JE)«UHI-U(JE) 

"  I F  (  DELUTJE  )  Vu  ( JE )  -UHT)  'Z~61*  266 . 206 
207  W(JE)»UL0-U(JE) 

_ I F  ( J_EJ+U (  JEJ.ry'-O  )  ^j6  *_2 06 j 208 
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208  W(JE)«DELU( JE) 

206  CONTINUE 

ETAN«0_.. 

DO  209  JE»1*KA 

209  ETAN»ETAN-W( JE)*C( JE) 

DO  211  I«1«N _ 

Z«0. 

DO  210  JE»1*ICA 

210  Z’Z+HII  »JEJ*W(JE) _  f 

211  DM(n*Z 
ETAD=0. 

_ PQ.  213-1  _ 

Z=0. 

DO  212  I«1.N 

212  Z«2-fP(12»I)*DM(l  ) _ 

1  =  12 

213  ETAD=ETAD+Z*DM( I ) 

_ EIA?i£lAN/ET.Afi _ j _ 

IF(1.“ETA)  214.216.216 

214  ETA*1. 

216  DO  2Jl 7  JE« Ij K A _ 

217  U(jn»U(JE)+'ETA*W(JE) 

250  ITER-l+ITER 

73  IF(MAX-ITER)971»971«77 _ 

77  DO  97  JE«1#ICA 

IF(C(JE))  81t802»82 

81  ABSC( JE)  =  -C(JE) _ 

GO  TO  802 

82  ABSC( JE)«C( JE) 

802  IF(U(  JE)-UHI  )  8^8^f84 _ 

84  WRITE  OUtPUT  T'aPE  3»85 

85  FORMAT (20H  UH  OUTSIDE  BOUNDARY) 

86  IF(ULO-U(  JE)  )  91»91»88 _ _ 

88  WRITE  OUTPUT  TAPE  3.89 

89  F0RMAT(20H  UL  OUTSIDE  BOUNDARY) 

91  IF(ABSC( JE)-EPSIL)  97.97.93 _ 

93  IF(C(JE))  94.97.95  , 

94  IF(UHI-U( JE)-DELTA)  97.97.96 

95  IF(ULO-U{ JE)+DELTA)  96.97.97 

96  GO'  to  20 

97  CONTINUE 
GO  TO  98 

971  WRITE  OUTPUT  TAPE  3*972 

972  FORMAT (21H  ITERATIONS  EQUAL  MAX) 

98  WRITE  OUTPUT  TAPE  3.99 
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99  FORMAT  {  35h”1  TERATl  ONs'  cOMP^  FOLLOWINg  Is^'u 

100  DO  101  JE  =  1*ICA 

WRITE  OUTPUT  TAPE  14,7*U(JE) 

101  WRITE  OUTPUT  TAPE  3»7.U(JE) 

WRITE  OUTPUT  TAPE  3.5. ITER 

_ 5_F0i^A1(  I7_) _  _ _  _ _ 

w’RITE  output  tape  3.7.V  . . 

CALL  EXIT 

_ END _ _ _  1 _ _  _ 
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The  firet  READ  etatement  aecapta  the  a;irier  of  the  systaat  the  nuibec  of 

control  aubintaxvalet  the  ■axlaieil  nunbar  of  IterationBf  the  upper  boundary 

for  control,  the  lower  boundary  for  control, and  the  values  for  £  and  P • 

Thaaa  are  N,  KA,  MAX,  UHI,  ULO,  ffSIL,  and  DELTA  reapectivaly.  KSENS 

Is  read  in  by  tha  next  atatamant.  A  KSENS  of  0  eelacta  tha  gradiant 

projection  daacent  achsaa  while  a  KSENS  of  1  aalacta  tha  "bang^ang" 

♦ 

daacent  schana*  p  is  next  read  in  by  atataaent  2  aa  A*  Statsaant  13 
aecepta  U  which  is  the  initial  gusas  for  the  control.  The  aatricaa  H 
and  R  are  read  in  by  stateaents  16  and  18.  Tha  DO  loop  ending  with 
statsaant  24  calculates  terainal  norm  V  is  calculated  by 

the  00  loop  ending  with  atataaent  35.  The  gradient  ia  calculated  next 
by  tha  DO  loop  with  atataaant  41.  Tha  next  atataaant  aalacta  tha  option 
of  using  the  Qradiant  Projection  Scheae  Or  the  Cbmtr  Aiming- Schaaa. 
Consider  the  Cemor-  AiMllng-^Sohaae  -firet*  -Tha'-DO  loop  beginning -with' fcf ate- 
asnt  42  and  ending  with  46  exaainea  each  coaponent  of  the  gradient 
vector  and  dstszaines  the  proper  according  to  thg  Comer  Aiming. Oesaint 
Scheae.  Statsaant  46  then  perfoxaa  the  aatrix  vector  aultiplication 
of  The  00  loops  ending  with  stateewnt  58  then  fora  the  product 

Tbe  next  statsaant  aalculatea  the  value. which  is  alwaya 
positive.  The  next  four  atateasnta  detezaine  Statsaant  72  than 

calculatsa  the  new  value  of  U.  Control  is  then  transferred  to  statsaant 
250. 

If  the  Sradiant  T>xcJ action  Descent  Scheae  had  been  selsctsd  control 
would  have  been,  taaaidbrred  from  statsaant  41'  to  statement  199.  The  DO 
loops  beginning  with  200  and  ending  with  the  atataaent  after  205  cal- 


The  00  loop  beginning  with  statement  S3  then 
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culate  '*1  *  ^  , 


selects  the  prdper  values  for  accordance  with  the  (gradient  Pro¬ 

jection  Scheme*  A  new  value  for  is  then  calculated  by  the  next 
sixteen  statements*  Statement  214  and  tha  one  preeadin^  it  ccMipute  /OtStOl)* 
The  00  loop  ending  with  217  then  conputse  the  new  value  for  U.  Statement 
250  adds  1  to  ITER,  which  is  the  number  of  iterations  completed.  The 
value  for  ITER  is  then  tested  to  datexmins  if  it  equals  MAX*  If  it 
does  the  value  of  U  Just  cosiputed  is  read  out*  If  ITER  dose  not  equal 
MAX  the  00  loop  beginning  with  statsmant  77  is  entered*  This  DO  loop 
detezmines  if  the  conditions  of  equation  (144)  for  the  termination 
of  the  descent  have  been  met*  If  not,  control  is  transferred  to  statement 
20  and  another  iteration  begins*  If  the  conditions  of  equation  (144) 
have  been  met,  the  next  few  statements  read  out  U,  ITER,  and  V*  The 
values  for  U  are  also  punched  on  cards  for  later  use* 


Calculation  of  V 

The  final  program,  which  is  labeled  CALCULATION  OF  V  calculates 
the  trajectories  of  the  state  variables*  This  program  which  is  listed 
on  the  next  two  pages  is  written  in  7090  FORTRAN,  It  requires  the 
matrix,  the  H  matrix  and  the  matrix  for  a  particular  problem!  it 
then  will  accept  as  many  different  solutions  for  the  control  ^  as 
are  desired  and  from  these  compute  the  trajectories  for  the  the  state 
variables  for  each  ^provided*  The  trajectories  are  calculated  by  solving 
for  the  values  of 

'KCt)  ^  ^  i-  e(^Z)  (158) 
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C  CALCOLATION  t)F  V - - - - - 

WRITE  OUTPUT  TAPE  3,20 

20  FORMATdSH  FOLLOWING' tS  V) 

READ  INPUT  TAPE  2,99,N,KA,KC  . .  ‘ 

99  FORMAT (31 6) 

DIMENSION  XT(6*eO) ,H(6*83) *U ( 80 ) ,p ( 6 *6 ) ,A( 6 *80 ) 

'  DO -I’  ME»1,»ICC  - - - 

DO  1  I«1,N 

1  read  input  tape  2, 18, A( I, ME) 

DO  2  JE*1,<A  ■  .  . ^ - 

DO  2  I«1,N 

2  READ  INPUT  TAPE  2,18,H(I,JE) 

. DO  6  '  - - -  - - - 

DO  6  I2«1,N 

6  READ  INPUT  TAPE  2,97,P(I2,I) 

97  FORMAT ( no. 0)  - - - 

31  READ  INPUT  TAPE  2,98,NB,KSENS 
K«1 

-'M-o'”  '  :  ■  : - 

DO  3  JB-1,NB 

3  READ  INPUT  TAPE  2,18,U{JB) 

'  9  ME«M+1  . . 

K8-NB-M 

IF(<C-ME)51,52,52 

'51"MEiKC  - - - 

52  DO  5  I«1,N 

2*0, 

. WA  - - - - — 

JE*JB+M+KA“NB 

4  Z*Z+H( I  * JE»*U( JB) 

'  -r-xTTHKr«z+A[Tvwrr - 

IF(l-K)  71,72,72 
72.  V*0. 

DO  -7' T2*r,N  - 

Z»0.  ' 

DO  8  I-1,N 

8  'z*z+xT(T  ;  iT*prrr»n - 

1*12 

7  V=V+Z*XT(  1,1) 

71  K*1+IC . -  - - - - 

M»M+1 

13  IF(NB-M) 11,11,9 

11  -WR  ITE  OOTPUT  TAPE  3 ','19  fV - 

41  WRITE  OUTPUT  TAPE  3,21 

21  FORMAT(30H  FOLLOWING  IS  XT  STARTING  AT't) 
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DO  16  fC  =  l»NB  . . . . 

16  WRITE  OUTPUT  TAPE  3 . 19  * ( XT ( 1 »<> » I «1 ♦ N ) 
IF(KSENS)31»31»30 
15  FORKAT(E10.4) 

18  FORMAT(E16.8) 

19  F0RMAT<6E12*4) 

98  FORMAT  (  1 8  *  14)  . . . .  . . 

30  CALL  EXIT 
END 
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at  each  control  subinterval*  In  order  to  compute  these  trajectories 
a  rather  involved  system  of  indexing  was  necassary. 

The  first  HEAD  statement  of  the  program  requires  the  order  of 
the  system,  the  number  of  columns  of  ^  and  the  number  of  columns  of 
In  case  constant  throughout  the  entire  control  interval  it  is 

only  necessary  to  supply  one  column  of  for  this  reason,  the  program 

was  designed  to  use  a  different  number  of  columns  of  f ^fjthan  of  /y  • 

The  READ  statements  down  to  statement  6  accept  the  matricas  H  » 

and  A.  .  Statement  31  requires  a  value  for  the  number  of  cards  to 
be  read  by  statement  3,  which  accepts  the  matrix  A  1  or  a  0  is 

also  required  by  statement  31  for  KSENS*  The  next  ten  statements 
down  to  statement  5  calculate  X(t)^Bach  element  of  which  is  a  point 
in  the  /trajectory  of  a  state  variable*  The  next  statement  tests  to  see 
whether  or  not  for  the  entire  control  interval  was  just  computed; 
if  it  was,  Vis  calculated  by  the  DO  loop  ending  with  statement  7; 

]((t)  was  not  just  calculated,  control  is  trsnsferred  to  statement  71* 

In  either  casa,  however,  new  indices  are  determined  by  71  and  the 
following  statement.  It  should  be  pointed  out  that  the  trajectories' 
are  calculated  beginning  atit*T  and  working  back  to  »  This  was  done 

to  prevent  further  indexing  complications*  Statement  13  tests  to  see 
if  all  the  points  in  the  trajectory  have  been  computed*  If  they  have 
not,  control  is  transferred  to  statement  9  for  further  computation* 

If  the  trajectories  are  complete  they  are  read  out  by  statement  16* 

At  this  point  the  value  for  KSENS  is  tested.  If  it  was  zero  there  will 
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be  another  set  of  trajectories  to  be  calculated  for  a  different  control 
interval;  thus,  control  will  be  transferred  to  31  where  new  values 
for  the  number  of  control  subintervals  and  for  KSENS  will  be  read  in. 

The  new  control  will  also  be  read  in  by  statement  3  and  the  computations 
then  repeated*  If  the  last  value  of  KSENS  was  1,  however,  exit  will 
be  called. 

Use  of  the  Programs 

The  use  of  the  programs  is  straightforeward  and  should  Jsresent 

0 

no  difficulty  if  the  following  points  are  kept  in  mind: 

1) '  all  matrices  are  read  in  by  columns;  each  element  of  a  matrix 
must  be  on  a  separate  card; 

2)  FORMAT  for  the  input  statements,  if  specified,  must  be  followed; 

3)  all  data  to  be  used  in  another ^program  is  punched  in  the  proper 
FORMAT  and  in  the  proper  order; 

4)  DIMENSION  statements  may  have  to  be  changed. 

Finally,  the  programs  may  be  modified  for  use  on  other  computers* 
However,  since  the  descent  program  contains  a  large  number  of  statements 
memory  limitations  of  the  computer  should  be  kept  in  mind. 
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Appendix  B 
Sample  Problem  #1 


In  order  to  check  out  the  computer  programs  a  sample  problem 
which  could  be  worked  out  by  hand  was  necessary.  It  was  felt  that 
if  an  interested  reader  wished  to  duplicate  the  results  of  this 
report  a  statement  of  this  example  problem  and  its  solution  might 
be  helpful.  Following  through  the  solution  will  also  give  a  better 
understanding  of  the  calculations  of  the  solutions  for  a  more  complex 
problem. 

The  problem  which  was  selected  is  a  very  simple  (from  the 
standpoint  of  calculation)  fourth  order  problem.  The  following  are 
the  system  equations  in  Stats  Space  form: 


^1  - 


^3 


•f- 


(159) 


A  control  interval  of  TtjjT  and  a  subinterval  ofZ*^  were  chosen. 
The  ^  and  H  matrices  can  be  calculated  by  direct  integration  of 
the  system  equations  as  was  developed  in  Chapter  IV.  To  find  the 
vector  (t,  the  system  equations  may  be  integrated  with  an  initial 
condition  of  The  system  must  be  excited  by  a  unit  step 
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input.  Thus,  to  find  jg,  integrate  from  tL:  O  to  ? the  system 


) 

l^,Lo)  -0 

TCy  =  -2i^  A  1 
^  TT  * 

^ 

~  ) 

y^CeO  ~o 

^  =  '*3  ; 

y-fCo)  =(^ 

(160) 


The  integration  yields  for  ^ 

z  C^;t4l  r  (161) 

^ 


Using  the  above  result  as  an  initial  condition  the  following  system 


can  be  integrated  from  zero  to  each  control  subinterval  in  order 
to  calculate  the  matrix ■  H  * 


■?</  s 

^ ; 

x,Co) 

- 

TCzC^') 

r  hfT- 

^3  '' 

=  .^93 

:r  .  107 

The  resuJ 

t  of  this  integration 

is  H  t 

•  tiy 

.0^/ 

■3fS 

~  .707 

•  Its 

.0/^8 

.  701 

-.39Z 

Its 

•0737 

.  707 

.  093 

.7SrS 

P<f3S' 

•393 

.  707 

.US 

.  ;a/ 

-.afs 

.  7  07 

H'- 

.78S 

./// 

~  .  707 

.3  73 

t  \ 

.78/ 

.19$ 

-.707 

^93 

.78/ 

.  35/ 

-•  39T . 

•'.707 

■7lS 

•  338 

.3f^ 

-  .  701 

.US 

.  930 

.  107 

-  •  2^5 

.78/ 

.S9^ 

■70? 

•393 

.71/ 

.  (,7f 

.  3  52 

•  709 

(162) 


(163) 
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Since  this  problem  is  used  to  check  results  of  the  computation  of 
the  optimum  driving  function,  an  optimum  driving  function  i^  chosen 
using  the  maximum  principal  and  then  the  initial  error  calculated 


by  using  this  driving  function,  and  sotting  the  solution  to  the 


This  correspond  to  a  A/u  of 

Q/  f  I  /  i]  (163) 

Forming  the  product  gives  the  vector 

\-sio 

Hx?'  -  y  (164) 

O  . 

The  matrix^^li^ust  now  be  calculated.  Again  the  system  equations 
may  be  integrated  directly  by  hand  to  compute  M  .  The  following 
is  the  matrix  equation  which  is  integrated  after  first  putting 
it  into  scalar  notation.  The  matrix  F  is  given  Table  B-It 

6£  ^  P  $  ;  =  I  (165) 

^  ±  - 
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The  following  is  obtained  for  ^  (  371 ) 

^(3Jl)  = 


10 

0 

O 

O 

0 

■OS' 

o 

O 

0 

0 

-1.0 

O 

0 

o 

D 

~l.t> 

(166) 


Now  set  the  solution  to  the  system  equal  to  zero 

£■+  H  ^  -  O 


Numerically  this  is 


t.oo  0  0 

s' 

jr 

0  .OS  o  o 

1.  S?o 

O  O  -t  o 

0  O  0  -l.O 

S-K 

...o 

The  initial  condition  C  is  readily  determined  to  be 

^JT 

^  30.^1  f 

i> 

O 


(167) 


(168) 


(169) 


Thus,  by  working  backward  from  a  known  solution  a  problem  has  been 
formulated  v;hich  can  be  used  to  test  the  digital  computer  programs. 
The  data  of  Table  B-I  was  used  for  the  computer  solutions. 

The  pertinent  results  of  the  computations  are  shown  in  Tables 
B-I I  and  B-I I I  on  the  last  page  of  this  appendix.  A  slight  change 
had  to  be  made  in  the  control  interval  and  subinterval  in  order  to 
make  subinterval  divide  into  the  control' interval  an  integral  number 
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r 

of  times.  This  change  accounts  for  the  slight  discrepancy  in  the 
third  decimal  place  of  the  ^  matrix.  The  descent  was  terminated 
after  135  iterations  when  it  was  evident  that  the  driving  function 
would  converge  to  that  calculated  above. 

Another  conclusion  may  be  drawn  from  the  results  of  this  problem 
in  addition  to  the  fact  that  the  computer  programs  were  correct. 

The  same  "bang-bang"  solution  whiph  was  used  to  determine  the  initial 
condition  was  arrived  at  by  an  entirely  ‘different  method  when  the 
computer  programs  were  used  to  find  the  driving  function.  The  maximum 
principal  was  used  to  formulate  the  control  in  the  first  case.  The 
optimization  techni'que  of  Steepest  Descent  was  used  in  the  second 
case.  This  result  again  supports  the  maximum  principal  of  Pontriagin. 


0 


o 


o 

-.319 

o 

0 

o 

0 

0 

1  ■ 

a. 

0 

'/ 

i2_ 

1 ' 

I 

-30.^11 

0 

(>.o 

J  _ 

O.o 

f 

0 

0 

0 

c> 

5 

o 

0 

0 

0 

2. 

0 

g_ 

0 

o 

/_ 

T 

Tfit/  -  -796 

uHf  = 

VLp  i  -/o 
£■  P$IL  -  -DOi 
QBLTA  -  .^Af 
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Table  B-II 


H  ^  Computed  by  Riinge-Kutta 


F0R  H  SW  0NE  0FF  ENTER  ALIT-AS  A 
F0R  LG  A  SW  0NE  0N  WNTER  XN0T  AS  A 
F0LL0WING  IS  H-TRANSP0SE 
'.7B6000E+00  '.V44853E-01 

'.7B6000E+00  -.57tl73E-01 

'.786000E+00 


-.786000E+00 

;7B6oooe+oo 
; 786000 E+00 
786000  E+00 
'.7860'00E+00 
'.7B6000E+00 
-.786000E+00. 
'.786000E+00 
.786000E+00 
T0P 


.;733363E-0T 
■.941607E-01 
'.1208S8E+00 
-.155228E+00 
■.199306E+00 
-.255901  E+00 
-.328567E+00 
-.421866E+00 
-.541658E+00 
.695467E+00 


-.288202  E+00 
;705754E+00 
-.709283  E+00 
-.29672OE+OO 
--.289909E+00 
-'.70  6468  E+00 
--.708586E+00 
-.295020 E+00 

-.291615E+00 

-.707I7BE+OO 

-.707884E+00 

.2933I8E+OO 


Table  B-II I 


Computed  Briving  Function 


1TER  EOUALS  135 
-.99IB78O6E+OO 
'.991B7B06E-+00 
'.99IB7B06E+OO  . 
'.99IB7B06E+OO 
-49866I603F-+OO 
-'.99I876O6E+OO 
-'.991 8780  6E+00 
--.9  63  9  8862  E+00 
-.991 8780  6E+00 
'.99I87806E+OU 
'.99I878O6F.+OO 
.99IB7806E+OO 


I 


-'.709631  E+00 
-'.29757OE+OO 
-..289056E+OO 
-.706112E+00 
-.70  893  5  E+00 
.295870 E+00 
-'.290762E+00 
-'.706824E+00 
-'.70  823  5  E+00 
-'.2941 69  E+00 
-.292467E+00 
.707531 E+00 
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Appendix  C 

Computer  Solution  of  a  Fourth  Order  Test  Problem 
In  order  to  further  test  the  computer  programs-  and  to  gain  . 
experience  in  their  use  a  more  realistic  problem  than  the  one  presented 
in  Appendix  B  was  used.  The  example  problem  solved  by  Ho  and  Brentani 
in  their  paper  (Ref  5:A-37)  was  chosen,  because  an  answer  was  readily 
availab.le.  Thtf  results  are  of  interest  in  that  a  "bang-bang"  solution 
was  again  obtained. 

The  data  for  the  programs  is  set  forth  in  Table  C-I.The  solution 
of  the  optimum  driving  function  for  a  control  interval  of  2*5  seconds 
is  presented  in  Table  These  results  agree  with  those  obtained 

by  Ho  and  Brentani*  ^ 
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G[t) 

F  r 

K  r 


Table  C-I 

Data  for  Computer  Programs 


(5+l)(Sf  a.)  C5*  +  f-  ^  ) 

0 

1 

0 

cT 

O 

0 

( 

0 

O 

0 

0 

7 

d- 

-10 

-10 

-£ 

~o 

.  5 

.  J  = 

P 

0 

l  A  V  •  OS 

1. 

C  : 

o 

l< 

-is 

o 

VLO  -  -/-O 

Bf^lL  '  oos. 

•  QBC.r/1=  025 

T 

0 

0 

0 

0 

3 

0 

o 

0 

0 

s. 

o 

o 

0 

0 

y  J 
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Table  C-II 

Optimum  Driving  Function 


ITER  EQUALS  50 
-.1 0000000 E+01 
-.tOOOOOOOE+01 
-.lOOOOOOOE+01 
-.10000000E+01 
-.lOOOOOOOE+01 
-.10000000E+01 
-.lOOOOOOOE+01 
-.1 0000000 E+01 
-.10000000E+01 
-,10000000E+01  . 
-.lOOOOOOOE+01 
-.10000000E+01 
•-.10000000E+01 
-.lOOOOOOOE+01  . 
-.lOOOOOOOE+01 
-..10000000E+01 
-.10000000E+01 
-,10000000E+01 
-.10000000E+01 
-.lOOOOOOOE+01 
-.100ODO0OE+O1 
-.lOOOOOOOE+01 
-.10000000E+01 

-.iooo6oooe+oi 

-JOOOOOOOE+01 

-.lOOOOOOOE+01 

-.lOOOOOOOE+Ol 

-.ionoooooE+01 

-.10000000E+01 

-.loooooooE+or 

-.10000000E+01 

-.lOOOOOOOE+Ol- 

-.lOOOOOOOE+01- 

-.loooooooE+or 

-.lOOOOOOOE+Ol' 

-•.3  1  93  2941E+00 
.69166105E+00 
.lOOOOOOOE+01 
.lOOOOOOOE+OT; 
.loooooooE+oi 

.lOOOOOOOE+01 
.lOOOOOOOE+Ol 
■  .lOOOOOOOE+Ol, 
.lOOOOOOOE+Ol! 

.  .loooobooE+or  • 

.35212415E+00, 

-.loooooooE+oi: 

r*.l0000000E+01 

-JOOOOOOOE+01 

-,10000000E+01 
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Appendix  D 
Tabular  DatA 

Tabular  data  from  the  computations  descrived  in  Chapter  V  is 
presented  in  this  appendix*  The  first  two  tables  are  the  matrices 
for  the  fourth  order  approximation  and  the  sixth  order  system  res¬ 
pectively.  The  rest  of  the  tables  are  optimum  driving  functions 
for  various  solution  and ' various  values  of  the  control  interval  7  • 
In  all  the  tables  of  the  optimum  driving  functions  the  last  two 
numbers  from  the  bottom  are  the  number  of  iterations  used  to  :oompute 
the  driving  function  shown  and  the  value  of  the  terminal  norm  V , 
respectively. 


130 


GE/EE/62-8 


'"FI)LL0WING  IS 

.ti?73yE-ot 
, .127419E-01 
..t3  533  6E-Ot 
.143445E-01 
i .t51695E-01 
.160031E-01 
.t68392E-01 
.17671.7E-01 
.t84938E-01 
.192985E-01  ^ 
.200784E-0t 
. 208261 E-01 
.2t  533  7E-Ot 
i .221931E-01 
. 227963 E-01 
.233351E-01 
.23  8013E-01 
.241867E-01 
. 24483 4E-Ot 
. 24683 4E-01 
.247793 E-01 
.24763 7 E-01 
.246298E-01 
.243714E-01 
.239825E-01 
.23 458 2 E-01 
.227940E-01 
.219864E-01 
.210326E-01 
.199311 E-01 
.186810E-01 
.172829E-01 
.157384E-01 
.140503 E-01 
.122228E-01 
. 10261 5E-01 
.817317E-02 
.59661 lE-0 2 
.385010E-02 
.123636E-02 


TablB  0-1 


Fourth  Oxdor  Approximation 


H  TRANSPOSE 

^.150925E-01 
-.156121E-01 
-.160431E^01 
-.163767E-01 
-.166043  E-01 
-.167174E-01 
-.167079E-01 
-.165684E-01 
-.162918E-Gt 
-.158716E-01 
-.153022E-0t 
-.145786E-01 
-.t36968E-ei 
-.12653  7  E-01 
-.tl4475E-Ot 
-.100773  E-01 
-.85^46E-02 
-. 684763  E-02 
-.499290E-02 
-.298374E-02 
-.826U4E-03 
.147 2 50 E-02 
.390307E-02 
. 64549 5E-02 
.911590E^02 
.1187 20 E-01 
.147079E-01 
.176065E-01 
.205492E-01 
.23  51  59E-01 
.264851 E-01 
.2943  3  7  E-01 
.3  23376E-01 
.351715E-01 
.379O89E-OI 
.405225E-01 
.429843 E-01 
.452655E-01 

.473373 E-01 
.49; 703 E-01 


.1121b2E-01 

.953558E-O2 

.767584E-02 

.564029E-02 

.343375E-O2 

.106267E-02 
-.t46484E-02 
-.413893  e-02 
-.694799E-02 
-.987866E-02 

-.'1291 57  E-01 

-.160423 E-01 
-.m397£-0t 
-.224876E-01 
-. 25763 8E-01 
-.290446 E-01 
-.323049E-01 
-.355181E-01 
-.386562E-Ot 
-.416902E-01 
-*445900E-01 
-.473 246 E-01 
-. 49862 5E-01 
-.521717E-OI 
-.542 200 E-01 
-.559749E-OI 
-.574047E-01 
-.584778E-01 
-. 59163 5E-01 
-.594324E-01 
-.592562E-01 
-.586084E-01 
-.574645E-01 
-.5580248-01 
-.53b026E-01 
-. 5084b 5E-01 
-.4752b7E-Ul 
-.436276E-G1 
-.391453  E-01 
-.340781  E-01 


.31C5091E-01 
.3541  25E-01 
.3b9662E-ul 
.424395E-U1 
.458OO5E-OI 
.490159E-01 
.52051fE-01 
.548731 E-01 
.57 4448 E-01 
.5973 11 E-01 
.616965E-01 
.633 056 E-01 
.64523 6 E-01 
.653  166  E-01 
.656518E-01 
.65498OE-OI 
.648258E-U1 
.63 6079 E-01 

.6l8l9eE-01 
.59439§E-U1 
.564486 E-01 
.5283 20 E-01 
.485786E-01 
.43  6  8  1  6E-01 
.381391  E-01 
.319533E-U1 
.251322E-01 
.I7689OE-UI 
. 964257 E-U2 
.101757E-U2 
-.815514E-02 
-.17838bE-01 
-.279896E-6I 
-.385586E-U1 
-.494o96£-01 
-.607204E-U1 
-.72182bE-01 
-.836OI6E-OI 
-.954970Et01 
-.107182E+Utj 
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Sixth  Order  Systeei 


S 
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m  m  n 
,000 
I  I  I 
ID  LU  LU 
r*-  ^ 

^  .CQ 
(\i  CC 
pO  •4’ 
fsj  rn  >r 


pnmf^fncgrgcgcsjfvi 

OOOOOOOOO 

t  I  I  I  (  I  I  I  I 

UJUJ'DDUJUJLUUJID 

r^>OCO^•^l^\G04'^• 

a)^co<“<Ni>t'00»-< 
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UJOIUJUJUJLUUJIDUJ 
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6E/EE/62-e 


Table  D>III 

OptiMum  Control  for  Sixth  Order  Syata« 
Direct  Control  of  Fiiot  Two  State  Variablaa 
T  ■  .675  aec* 


■  ITtRATIONS  COMPLETE  FOLLUWINu  IS  U 

-0.  IS»500000E 

03 

-0. 19 500000 E 

03 

-0.19500000E 

03 

-0.19500000E 

03 

-U.  19500000F 

03- 

'  -0.19500000E 

03 

-0. 195000()0E 

03 

-0.19500000E 

03 

-0.  l9500U00r 

03 

-0.19500000E 

03 

-0. 19500000b 

03 

-0.19500000E 

03 

-0. 19500U00E 

03 

-0.19500000E 

03 

-0.19500000E 

03 

-0.19500000E 

03 

-C.IU524182F 

03. 

0, 19500000C 

03 

0.19500000E 

03 

0. 19500000E 

03 

0. 19500000E 

03 

0.19500000E 

03 

0.19500000E 

03 

0.19500000E 

03 

p.l9500g00E 

03 

0.19500000E 

03 

0.19500000C 

03 

129 

0.  l203Qg42E- 

-00 
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GE/EC/62-8 


Table  D-III 
(eont.L 

Optimum  Control  for  Sixth  O^^er  S^tam 
Direct  Control  of  Pint  Two  State  Variablea 
T  ■  *625  eae* 


1  ITERATIONS  COMPLETE  FOLLOW INo  IS  U 

_-0.19500000E 

03 

-0.19500000F 

03 

-0. 19500000E 

03, 

-0.19500000E 

03 

-0.1950d000E 

03 

-0.19500000E 

03 

-0.19500000E 

03 

-0. 1950000CE 

03 

-0.19500UOOE 

03 

-0. 19500000E 

03 

-0.19500000E 

03  ■ 

-0. lybOOOOOt 

03 

-0.19500U00E 

03 

-0.l9iiOOOOOE 

03 

-0.19500000E 

03 

-0.12865712E. 

02 

0.19500000E 

03 

0.19500000E 

03 

Q.19500000E 

03 

0. 19500000E 

03 

0.1950000UE 

03 

0.19500()OOE 

03 

0.19500000E 

03  . 

0. 19500000E 

03 

0.195000QQE 

03 

60 

0.35551506E- 

-00  : 
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Ge/eE/62-6 

Table  0-  III 
(cent*) 

OptiaMi  Control  for  Sixth  Order  Syatan 
Oisaet  Control  of  Pixat  Two  Stila  Variablaa 
T  ■  .575  aoe. 


irERAflONS  COMPLETE  FOLLOWING  IS  U 

„rJl-.195a.QaQQE 

DJ 

-C. lySOOOOOE 

03 

>0.19500000E 

03 

-0. 19500000E 

03 

-0.19500000E 

03 

^0. 19500000E 

03 

.. -D.  1.93QQQQ0E 

03, 

-U. 19500000E 

03 

-0. 19500000E 

03' 

-U. 19500000E 

03 

-0.19500000.E 

03 

-0. 19500000E 

03 

.  .-e..J.35fiQQQ0u 

03 

-0. 19500000E 

03 

0.95030949E 

02 

0.19500000E 

03 

0.19500000E 

03 

.0.19500000E 

03 

..  o.a9.;>oooooE 

03 

0.19500000E 

03 

0,l95ODOOOE 

03 

L. 19500000 E 

03 

0.195OOUOOE 

03 

65 

0.69346071E 

00 
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6C/EC/62-6 


Tabl*  D-III 
(eont*) 

OptiJMjM  Control  for  Sixth  Ordor  Systaa  . 
Diroet  Control  of  First  Two  Stats  Variables 
T  ■  «S00  sac* 


iTCRATlUNS  COMPLfeTC  FOLLOWING  IS  U 

-G.19bG0C)00E 

03 

-0.19500000E 

03 

-G. 19500000E 

03 

19500000E 

03 

-0.19900000E 

03 

-0. 19500000E 

03 

-0. 195000QOE 

03 

-0.1950000dE 

03 

,-0. 19500000E 

03 

-0.19500000E 

03 

-0,19500000E 

03 

-0.19500000E 

03 

Q.5^AflAeL40E 

02 

0.>9500000E 

03 

_ £L.  1.95011000  E 

03 

t).l9b00000E 

03 

..  0.^19500000 E 

03 

0.19500000E 

03 

_ O.lOSOOOOflE. 

-03  . 

0. 19500000E 

03 

44 

0.13508035E. 

01 

6E/Ee/<2^ 


Table  D-III 
(cent* ) 

Optiaun  Control  for  Sixth  Ordar  Syatan 
Diraet  Control  of  first  Two  Ststa  Variables 
T  ■  .375  sac. 


ITERATIONS  CUHPLETe  FOLLOWING  IS  U 

-&.195OOO006 

03 

-0. 19500000E 

03 

-0.l91j00000E 

03 

-0.l9b00000E 

03 

-0-19500UO0E 

03 

-0.19500000E 

03 

-0,19300000E 

03 

-0. 19500000E 

03 

-0.19i>00000E 

03 

0.18410978E 

03 

0-19500000E 

03 

G-195()0000E 

03 

0.19500000E 

03 

0.19b00000E 

03 

0.l9i)OOOOOE 

03 

33 

0.26502431E 

01 

Se/EE/62^ 


Tabla  D>  IV 

OptlMM  Conitol  for  Fourth  Ordor  ApproxiMtian 
Dlxaet  Control  of  Fimt  Two  State  Variablaa 
T  >  .650  aac. 


ITeKATIQNS  COMPLETE  FOLLOWING  IS  U 
-0.7U999999E  01 
-0.70999999E  01 . 

-0,70999999E  01 
-0.70999999E  01 
-Q. 70999999E  01 
-0.70999999E  01 
.-p,70999999E.  01 
’-0.709q'9999E  01 
0.23082887E  01 
0.70999999E  01 
0.70999999E  01 
0.709q9999E  01 

0.70999999E*01  _  _ 

62  . . 

,0.16079205E;00 _ 


T  -  .600  aae. 


ITERATIONS  COMPLETE  FOLLOWING  IS  U 

^.Ifl999999E  OX  _ 

-0.70999999E  01 

-zfl.7Q999939£.pL  _ 

-0.70999999E  01 
-0.70999999E  01 
-0.70999999E  61 

-Q.70999999E  01  _ 

-U.35306581E  01 
■Q.7Q999999E  01 
0.70999999E  01 
Q.70999999E  01 
0.70999999E  01 

_ Z1  _ _ , _ _ 

0.A1684502E-00 
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&E/EE/62-a 


Tablfl  0-  IV 

OptiMM  Control  for  Fourth  Ordor  Approximation 
Dirset  Control  of  First  Two  Stata  Variablea 
T  ■  «50Q  sac. 


ITERATIONS  C0^'PLETE  FOLLGWING  IS  U 

-0.7C999999E 

01 

-0.70999999E 

Cl 

-0.70999999E 

Cl 

-0.70999999t 

Cl 

-0.70999999E 

01 

-0.70999999E 

lOl 

‘  0.82834515E 

CC 

0.7C999999E 

01 

0.70999999E 

Cl 

0.70999999E 

01 

32 

0.121S8052E 

01 
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CE/CE/62-8 


Tobl*  D- V 

OptiiM  Contxol  for  Fourth  Order  ApproxiMotiori 
Direct  Centrol  of  All  Four  Stete  Veridblee  ' 
T  ■  1.85  etc. 


i  llTERAltpNS  CCHPLETE  FOLLOWING  IS 
'  -b.709999<J9E  01 
-6.7C999999E  01 
-0.7C999999E  01 
-0.709999996  01  i 

-6.709999996  01 
-0.709999996  01 
-6.709999996  01 

-0;76999999e  01  ,  ' 

-6.70999999E  01 

-6. 709999996;  01 

-6. 70999.9996  01 

-6.7099999:96  01 

-0.7C999999E  01 


0:, 


-0.7C999999E  01 
-0.709999996  01 
-0.6S408815E  01 
0“.709999996  01 
6.709999996  01 
0,;^ 709999996  01 
0.709999996  01 
0.7C999999E  01 
0‘.7C9999996  01 
0.709999996  01 
0.709999996  Cl 
__g» 709999996 ’01 
0.709999996  01 
0. 709999996  .01 
0.709999996  01 
-0.3A1565816  01 
-0.709999996  01 
_-0i70999999E  01 
-6i7C999999E  01 
-0. 709999996  01 
-0.709999996  01 
-0. 709999996  01 
-0:70999999E  01 
-0.709999996  01 
1350 

0.816957236-01 
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GE/EE/62o8 


Table  D*V 
(cont.) 

]ptiJMM  Conttol  for  Fourth  Ordar  ApproxiMation 
Direct  Control  of  Ail  Four  State  Variablea 
T  ■  1.75  sec. 

ITEKATIONS  CUMPLFTE  FOLLOk I NG. I S  U 
~0.70999999F  0 1' 

-0.7C999999E  01 
-0.70999999E  01 
-0.7U999999E  01 
-0.70999999fc  01 
-0.7C999999F  01 
-().70999999E  01 
'0.7U999999E  01 
-0.70999999E  01 
-0.7G9g99g9F  01 
-0.70999999E  01 
-0. 7U999999E  01 
-U.7C999999E  01 
-0.7C999999E  01 
-0. ltf88l38lE  01 
0.70999999E  01 
0. 70999999E  01 
0.70999999E  01 
0.70999999E  01 
0.70999999E  01 
C.7C999999F  01 
0.70999999E  01 
0.70999999E  01 
0.70999999E  01 

..  Q.tl&9999gqf  gi 

0.70999999E  01 
0,340A'3817E  01 
-0.70999999E  01 
-0.70999999E  01 
-0.70999999E  01 
-Q.-IQ999.9.99E  01 
-0.70999999E  01 
-0.70999999E  0,1 
-0.70999999E  01 
-0.709999"99E  01 
69^> 

0.2^e7U2,5E-00 
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GC/CC/62«e 


Table  D>V 
(cent.) 

OptlM*  Contxol  for  Fourth  Order  Approximation 
Direct  Control  of  All  Four  State  Variablea 


T  m  1,50  sac. 


irtRATIONS  COMPLETE  FOLLOb^ING  IS  U 

-0.70999999E 

01 

-0.70999999E 

01 

-0.70999999E 

01 

-0. 70999999E 

01 

-0.70999999E 

01 

-0.70999999E 

01 

-0.70999999E 

01 

-0. 70999999E 

01 

-0.70999999E 

01 

-0.70999999E 

01 

-0,70999999E 

01 

-0.22782274E 

01 

0.70999999E 

01 

0.70999999E 

01 

6.70999999E 

01 

0.70999999E 

01 

0.70999999E 

01 

0.70999999E 

01 

0,70999999E 

Qi 

0.70999999E 

01 

0.70999999E 

01 

0.70999999E 

01 

0.67485256E 

01 

-0,70999999E 

01 

-0.70999999E 

01.. 

-0.70999999E 

01 

-0.7 09999 99E 

01 

-0.70999999E 

01 

-0.70999999E 

01 

-0.7()999999t 

01 

.  355 

0.90P23644E 

00 
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GE/EE/62-8 


Table  D«V 
(cont.) 

Optimum  Control  for  Fourth  Order  Approximation 
Direct  Control  of  All  Four  State  Variablee 


T  ■  1.25  8BC. 


ITERATIONS  COMPLETE  FOLLOWING  IS  U 

-0.70999999E 

01 

-0.70999999E 

01 

-0.70999999E 

01 

-0.70999999E 

01 

-0.70999999E 

01 

, -0.70999999E 

01 

-0.70999999E 

01 

-0.70999999E 

01 

-0.70999999E 

01 

0.55A72609E 

01 

0.70999999E 

01 

0.70999999C 

01 

0.70999999E 

01 

0.70999999E 

01 

0.70999999E 

01 

0.70999999E 

01 

0.70999999E 

01 

0.70999999E 

01 

0.70999999E 

01 

-0.39656283C 

01 

-0.70999999E 

01 

-0.70999999E 

01 

-0.70999999E 

01 

-0.70999999E 

01 

-0..70999999E 

01 

290 

0.16630284E 

01 

GE/EE/62-8 


Table  D*V 
(eont*) 

OptijMja  Control  for  Fourth  Order  Approximation 
Direct  Control  of  All  Four  State  Variablee 
T  ■  l.OQ  eec* 


ITERATIONS  COMPLETE  FOLLOWING  IS  U 
-0.709999196  01 
-0.70999999E  01 
-0.709999996  01 
-0.70999999E  01 
-0.70999999E  01 
-0.70999999E  01 
-0.6883Q354E  01 
0.70999999E  01 
0.709999996  01 
0.709999996  01 
0.709999996  01 
0.709999996  01 
0.70999999E  01 
0.70999999E  01 
0.709999996  01 
0.U569404E-00 
-0.709999996  01 
-0.709999996  01 
-0.709999996  01 
-0.709999996  01 
460 

0.233949926  01 


< 
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Table  D-V 
(cont.) 

Optimum  Control  for  Fourth  Order  Afiproximetion 
Direct  Control  of  All  Four  State  Variablee 


T  m  .75  aec* 


ITERATIONS  COMPLETE  FOLLOWING  IS  U 

-0.70999999t 

01 

-0.70999999E 

01 

-0.70999999E 

01 

-0.70999999E 

01 

-0.70999999E 

01 

0.36115248E 

01 

0.70999999E 

01 

0.70999999E 

01 

0.70999999E 

01 

0.70999999E 

01 

0.70999999E 

01 

0.53904514E 

01 

-0.70999999t 

01 

-0.70999999E 

01 

-0.70999999E 

01 

337 

6.28406490E 

01 

T  ■  .50  see. 


ITERATIONS  COMPLETE  FOLLOWING  IS  U 

-0.70999999E 

01 

-0.70999999E= 

01 

-0.70999999E 

01 

-0.52056129E 

01 

0.70999999E 

01 

0.70999999E 

01 

0.70999999E 

01 

0.70999999E 

01 

-0.19497150E 

01 

-0.70999999E 

01 

754 

0.31596032E 

01 
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T^la  D-VII 

Stata  Vaxlabla  Trajactorias 
Fourth  Ordar  /|pp£oxiMation 
Dixaet  Contret  of  Two  Stata  Varidblaa 


FOLLCWING  IS 

XT  STARTING 

AT  T 

0.1582E-CC 

-0.178CE- 

GO 

0.  1994E 

01 

-0.1 8806 

01 

0.17866-00 

-0. 64036 

00 

0.22856 

01 

-0.83686 

00 

O.-2224E-0C 

-C.1116E 

01 

0.2521E 

01 

C.3123E- 

00 

0.2<3036-CC 

-0.160  IE 

01 

0.2698E 

01 

0. 1566E 

01 

0.38276-00 

-6.2094E 

01 

0.28086 

01 

' 0.2920E 

01 

0.4938E-0C 

-0.23486 

ul 

0.2713E 

01 

0.38026 

01 

0.6C55E  00 

-C.2118E 

01 

0.2305E 

01 

C. 36036 

01 

O.7053E  00 

-0. 1869E 

Cl 

0. 19C96 

01 

C.3329E 

01 

0. 79216  00 

-0. 1599E 

01 

0.15286 

01 

0.2978E 

01 

0.86496  CO 

-0.13126 

01 

0. 1167E 

01 

0.2547E 

01 

0.9229E  CO 

-C.1CC7E 

01 

0.8296E 

Ob' 

. C.2C34E" 

01 

0.96536  00 

-0.6852E 

00. 

0.5199E 

00 

C.  1439E 

01 

0.99126  CO 

-0.34916- 

00 

0.2420E- 

00 

0.7610E 

00 
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&.vin 

St«te  Varii*!*  TMjwtoriaa 
foiiirth  Onlir  ApfraexiMtian 
OlxKt  Contxel  of  Four  State  Varitelaa 


FOLLOWING  IS  XT  SIAkTING 

AT  T 

0.1221E-00 

-0. 7  0  38 E 

-03 

0. b481E-01 

-0.8  755E- 

-01 

0.  U33C-00 

0.35!52E- 

-00 

-0. 1 129E-00 

-0.9223E 

00 

0.8645E-01 

0-7199E 

00 

-0.2668E-00 

-0. 1837E 

01 

0.4120t-01 

0. 109ie 

01 

-0.3731E-00 

-0.2o24F 

01 

-0.2272e-01 

0. 1466E 

01 

-0.4278E-00 

-0.3896E 

01 

-0. 1055E-00 

0. 1843E 

01 

-0.4275E-00 

-0.8033E 

01 

-0.2070E-00 

0.2218E 

di 

-0.36fi6E-00 

-0.6238E 

01 

.-0.3272E-00 

Oi,2589E 

01 

-0.2478E-00 

-0.7‘304E 

01 

-0.4657E-Q0 

0.2952E 

01 

-0.6238E-()1 

-0.8826E 

01 

•-0.6175E  00 

0.3117E 

01 

0.2951E-00 

-0.9759E 

01 

-0.7638E  00 

0.2730E 

01 

b.9‘J83E.  00 

-0.9442E 

01 

-0.8899E  00 

0.2307E 

01 

0.1682E  01 

-0.8998E 

01 

-0.9940E  00 

0il852E 

01 

0.2341E  01 

-0.8421E 

01 

.-0.1075E  01 

0. 1364E 

01 

0.2967E  01 

-0.7708E 

01 

-O.inOE  01 

0.8456E 

00 

0.3I>55E  oi 

-0.6856E 

01 

“0.1159E  01 

e.2992E- 

-00 

0.4096E  01 

-O.5063E 

01 

-0.1i59E  01 

-0,2731E 

-00 

6.4f>84E  01 

-0.4727E 

01 

.0.1 

-0.1072E  01 

-0.1483E 

01 

0.5371E  01 

-0.2030E 

01 

-0.9824E  00 

-0.J114E 

01 

0.5657E  01 

-0.4713E- 

-00 

-O.06O6E  00 

-0.2^J7E 

01 

0.5b62E  oi 

6.i22'3E 

01 

.r!,G.7065E  OQ 

-0.3409E 

01 

0.5970E  01 

0.3049E 

01 

r0.5369E  00 

-0.3372E 

oi 

0.5612E  01 

0.3379E 

01 

j=Slx21Q^i-QQ. 

^Qj,.3287E 

01 

.e..5J15„E  01. 

0,.  3(i06E. 

QJ 

-0.2085E-00 

-0.3183E 

oi 

0.4807E  01 

0.3785E 

01 

-0.5241E-01 

^0. 3058E 

01 

0.4i92E  01 

0.3910'E 

01 

6.9692E-01 

-0.2912E 

01 

0.3972E  01 

0.3979F 

01 

0.2384E-00 

-0.2745E 

01 

O.3550E  01 

0.3986E 

01 

0.3711E-00 

-0.2557E 

01 

0.3129E  01 

0.3928E 

01 

-.Qt^.93^r00 

01 

0.2.713E  01 

0.3802F 

01 

0.6C55E  00 

-0.2118E 

01 

0.2305E  01 

d.36d3E 

01 

0.7053E  00 

^0.1869E 

01 

0.1909E  01 

0.3329E 

01 

0.7921E  00 

-0. 1599E 

01 

d.lb28E  01 

0.2978E 

01 

0.8649E  00 

-0. 1312E 

01 

0.1167E  01 

0.2547E 

01 

0.9229E  00 

-0. 1007E 

01 

0.8296E  00 

0.2034E 

01 

0.9653E  00 

-0.6852E 

00 

0,5199E  00 

0.1439E 

01 

.,.0.9912E  00 

-0.3491E- 

-00 

0.2420E-00 

■0.761PE 

00 

148 


